1 #include "blas_extended.h"
2 #include "blas_extended_private.h"
BLAS_cgemv2_c_s_x(enum blas_order_type order,enum blas_trans_type trans,int m,int n,const void * alpha,const void * a,int lda,const float * head_x,const float * tail_x,int incx,const void * beta,void * y,int incy,enum blas_prec_type prec)3 void BLAS_cgemv2_c_s_x(enum blas_order_type order, enum blas_trans_type trans,
4 int m, int n, const void *alpha, const void *a,
5 int lda, const float *head_x, const float *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 void*
34 *
35 * lda (input) int
36 * Leading dimension of A
37 *
38 * head_x
39 * tail_x (input) const float*
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_cgemv2_c_s_x";
62 switch (prec) {
63 case blas_prec_single:{
64
65 int i, j;
66 int iy, jx, kx, ky;
67 int lenx, leny;
68 int ai, aij;
69 int incai, incaij;
70
71 const float *a_i = (float *) a;
72 const float *head_x_i = head_x;
73 const float *tail_x_i = tail_x;
74 float *y_i = (float *) y;
75 float *alpha_i = (float *) alpha;
76 float *beta_i = (float *) beta;
77 float a_elem[2];
78 float x_elem;
79 float y_elem[2];
80 float prod[2];
81 float sum[2];
82 float sum2[2];
83 float tmp1[2];
84 float tmp2[2];
85
86
87 /* all error calls */
88 if (m < 0)
89 BLAS_error(routine_name, -3, m, 0);
90 else if (n <= 0)
91 BLAS_error(routine_name, -4, n, 0);
92 else if (incx == 0)
93 BLAS_error(routine_name, -10, incx, 0);
94 else if (incy == 0)
95 BLAS_error(routine_name, -13, incy, 0);
96
97 if ((order == blas_rowmajor) && (trans == blas_no_trans)) {
98 lenx = n;
99 leny = m;
100 incai = lda;
101 incaij = 1;
102 } else if ((order == blas_rowmajor) && (trans != blas_no_trans)) {
103 lenx = m;
104 leny = n;
105 incai = 1;
106 incaij = lda;
107 } else if ((order == blas_colmajor) && (trans == blas_no_trans)) {
108 lenx = n;
109 leny = m;
110 incai = 1;
111 incaij = lda;
112 } else { /* colmajor and blas_trans */
113 lenx = m;
114 leny = n;
115 incai = lda;
116 incaij = 1;
117 }
118
119 if (lda < leny)
120 BLAS_error(routine_name, -7, lda, NULL);
121
122
123
124
125 incy *= 2;
126 incai *= 2;
127 incaij *= 2;
128
129 if (incx > 0)
130 kx = 0;
131 else
132 kx = (1 - lenx) * incx;
133 if (incy > 0)
134 ky = 0;
135 else
136 ky = (1 - leny) * incy;
137
138 /* No extra-precision needed for alpha = 0 */
139 if (alpha_i[0] == 0.0 && alpha_i[1] == 0.0) {
140 if (beta_i[0] == 0.0 && beta_i[1] == 0.0) {
141 iy = ky;
142 for (i = 0; i < leny; i++) {
143 y_i[iy] = 0.0;
144 y_i[iy + 1] = 0.0;
145 iy += incy;
146 }
147 } else if (!(beta_i[0] == 0.0 && beta_i[1] == 0.0)) {
148 iy = ky;
149 for (i = 0; i < leny; i++) {
150 y_elem[0] = y_i[iy];
151 y_elem[1] = y_i[iy + 1];
152 {
153 tmp1[0] = y_elem[0] * beta_i[0] - y_elem[1] * beta_i[1];
154 tmp1[1] = y_elem[0] * beta_i[1] + y_elem[1] * beta_i[0];
155 }
156
157 y_i[iy] = tmp1[0];
158 y_i[iy + 1] = tmp1[1];
159 iy += incy;
160 }
161 }
162 } else { /* alpha != 0 */
163 if (trans == blas_conj_trans) {
164
165 /* if beta = 0, we can save m multiplies:
166 y = alpha*A*head_x + alpha*A*tail_x */
167 if (beta_i[0] == 0.0 && beta_i[1] == 0.0) {
168 if ((alpha_i[0] == 1.0 && alpha_i[1] == 0.0)) {
169 /* save m more multiplies if alpha = 1 */
170 ai = 0;
171 iy = ky;
172 for (i = 0; i < leny; i++) {
173 sum[0] = sum[1] = 0.0;
174 sum2[0] = sum2[1] = 0.0;
175 aij = ai;
176 jx = kx;
177 for (j = 0; j < lenx; j++) {
178 a_elem[0] = a_i[aij];
179 a_elem[1] = a_i[aij + 1];
180 a_elem[1] = -a_elem[1];
181 x_elem = head_x_i[jx];
182 {
183 prod[0] = a_elem[0] * x_elem;
184 prod[1] = a_elem[1] * x_elem;
185 }
186 sum[0] = sum[0] + prod[0];
187 sum[1] = sum[1] + prod[1];
188 x_elem = tail_x_i[jx];
189 {
190 prod[0] = a_elem[0] * x_elem;
191 prod[1] = a_elem[1] * x_elem;
192 }
193 sum2[0] = sum2[0] + prod[0];
194 sum2[1] = sum2[1] + prod[1];
195 aij += incaij;
196 jx += incx;
197 }
198 sum[0] = sum[0] + sum2[0];
199 sum[1] = sum[1] + sum2[1];
200 y_i[iy] = sum[0];
201 y_i[iy + 1] = sum[1];
202 ai += incai;
203 iy += incy;
204 } /* end for */
205 } else { /* alpha != 1 */
206 ai = 0;
207 iy = ky;
208 for (i = 0; i < leny; i++) {
209 sum[0] = sum[1] = 0.0;
210 sum2[0] = sum2[1] = 0.0;
211 aij = ai;
212 jx = kx;
213 for (j = 0; j < lenx; j++) {
214 a_elem[0] = a_i[aij];
215 a_elem[1] = a_i[aij + 1];
216 a_elem[1] = -a_elem[1];
217 x_elem = head_x_i[jx];
218 {
219 prod[0] = a_elem[0] * x_elem;
220 prod[1] = a_elem[1] * x_elem;
221 }
222 sum[0] = sum[0] + prod[0];
223 sum[1] = sum[1] + prod[1];
224 x_elem = tail_x_i[jx];
225 {
226 prod[0] = a_elem[0] * x_elem;
227 prod[1] = a_elem[1] * x_elem;
228 }
229 sum2[0] = sum2[0] + prod[0];
230 sum2[1] = sum2[1] + prod[1];
231 aij += incaij;
232 jx += incx;
233 }
234 {
235 tmp1[0] = sum[0] * alpha_i[0] - sum[1] * alpha_i[1];
236 tmp1[1] = sum[0] * alpha_i[1] + sum[1] * alpha_i[0];
237 }
238
239 {
240 tmp2[0] = sum2[0] * alpha_i[0] - sum2[1] * alpha_i[1];
241 tmp2[1] = sum2[0] * alpha_i[1] + sum2[1] * alpha_i[0];
242 }
243
244 tmp1[0] = tmp1[0] + tmp2[0];
245 tmp1[1] = tmp1[1] + tmp2[1];
246 y_i[iy] = tmp1[0];
247 y_i[iy + 1] = tmp1[1];
248 ai += incai;
249 iy += incy;
250 }
251 }
252 } else { /* beta != 0 */
253 if ((alpha_i[0] == 1.0 && alpha_i[1] == 0.0)) {
254 /* save m multiplies if alpha = 1 */
255 ai = 0;
256 iy = ky;
257 for (i = 0; i < leny; i++) {
258 sum[0] = sum[1] = 0.0;;
259 sum2[0] = sum2[1] = 0.0;;
260 aij = ai;
261 jx = kx;
262 for (j = 0; j < lenx; j++) {
263 a_elem[0] = a_i[aij];
264 a_elem[1] = a_i[aij + 1];
265 a_elem[1] = -a_elem[1];
266 x_elem = head_x_i[jx];
267 {
268 prod[0] = a_elem[0] * x_elem;
269 prod[1] = a_elem[1] * x_elem;
270 }
271 sum[0] = sum[0] + prod[0];
272 sum[1] = sum[1] + prod[1];
273 x_elem = tail_x_i[jx];
274 {
275 prod[0] = a_elem[0] * x_elem;
276 prod[1] = a_elem[1] * x_elem;
277 }
278 sum2[0] = sum2[0] + prod[0];
279 sum2[1] = sum2[1] + prod[1];
280 aij += incaij;
281 jx += incx;
282 }
283 sum[0] = sum[0] + sum2[0];
284 sum[1] = sum[1] + sum2[1];
285 y_elem[0] = y_i[iy];
286 y_elem[1] = y_i[iy + 1];
287 {
288 tmp1[0] = y_elem[0] * beta_i[0] - y_elem[1] * beta_i[1];
289 tmp1[1] = y_elem[0] * beta_i[1] + y_elem[1] * beta_i[0];
290 }
291
292 tmp2[0] = sum[0] + tmp1[0];
293 tmp2[1] = sum[1] + tmp1[1];
294 y_i[iy] = tmp2[0];
295 y_i[iy + 1] = tmp2[1];
296 ai += incai;
297 iy += incy;
298 }
299 } else { /* alpha != 1, the most general form:
300 y = alpha*A*head_x + alpha*A*tail_x + beta*y */
301 ai = 0;
302 iy = ky;
303 for (i = 0; i < leny; i++) {
304 sum[0] = sum[1] = 0.0;;
305 sum2[0] = sum2[1] = 0.0;;
306 aij = ai;
307 jx = kx;
308 for (j = 0; j < lenx; j++) {
309 a_elem[0] = a_i[aij];
310 a_elem[1] = a_i[aij + 1];
311 a_elem[1] = -a_elem[1];
312 x_elem = head_x_i[jx];
313 {
314 prod[0] = a_elem[0] * x_elem;
315 prod[1] = a_elem[1] * x_elem;
316 }
317 sum[0] = sum[0] + prod[0];
318 sum[1] = sum[1] + prod[1];
319 x_elem = tail_x_i[jx];
320 {
321 prod[0] = a_elem[0] * x_elem;
322 prod[1] = a_elem[1] * x_elem;
323 }
324 sum2[0] = sum2[0] + prod[0];
325 sum2[1] = sum2[1] + prod[1];
326 aij += incaij;
327 jx += incx;
328 }
329 {
330 tmp1[0] = sum[0] * alpha_i[0] - sum[1] * alpha_i[1];
331 tmp1[1] = sum[0] * alpha_i[1] + sum[1] * alpha_i[0];
332 }
333
334 {
335 tmp2[0] = sum2[0] * alpha_i[0] - sum2[1] * alpha_i[1];
336 tmp2[1] = sum2[0] * alpha_i[1] + sum2[1] * alpha_i[0];
337 }
338
339 tmp1[0] = tmp1[0] + tmp2[0];
340 tmp1[1] = tmp1[1] + tmp2[1];
341 y_elem[0] = y_i[iy];
342 y_elem[1] = y_i[iy + 1];
343 {
344 tmp2[0] = y_elem[0] * beta_i[0] - y_elem[1] * beta_i[1];
345 tmp2[1] = y_elem[0] * beta_i[1] + y_elem[1] * beta_i[0];
346 }
347
348 tmp1[0] = tmp1[0] + tmp2[0];
349 tmp1[1] = tmp1[1] + tmp2[1];
350 y_i[iy] = tmp1[0];
351 y_i[iy + 1] = tmp1[1];
352 ai += incai;
353 iy += incy;
354 }
355 }
356 }
357
358 } else {
359
360 /* if beta = 0, we can save m multiplies:
361 y = alpha*A*head_x + alpha*A*tail_x */
362 if (beta_i[0] == 0.0 && beta_i[1] == 0.0) {
363 if ((alpha_i[0] == 1.0 && alpha_i[1] == 0.0)) {
364 /* save m more multiplies if alpha = 1 */
365 ai = 0;
366 iy = ky;
367 for (i = 0; i < leny; i++) {
368 sum[0] = sum[1] = 0.0;
369 sum2[0] = sum2[1] = 0.0;
370 aij = ai;
371 jx = kx;
372 for (j = 0; j < lenx; j++) {
373 a_elem[0] = a_i[aij];
374 a_elem[1] = a_i[aij + 1];
375
376 x_elem = head_x_i[jx];
377 {
378 prod[0] = a_elem[0] * x_elem;
379 prod[1] = a_elem[1] * x_elem;
380 }
381 sum[0] = sum[0] + prod[0];
382 sum[1] = sum[1] + prod[1];
383 x_elem = tail_x_i[jx];
384 {
385 prod[0] = a_elem[0] * x_elem;
386 prod[1] = a_elem[1] * x_elem;
387 }
388 sum2[0] = sum2[0] + prod[0];
389 sum2[1] = sum2[1] + prod[1];
390 aij += incaij;
391 jx += incx;
392 }
393 sum[0] = sum[0] + sum2[0];
394 sum[1] = sum[1] + sum2[1];
395 y_i[iy] = sum[0];
396 y_i[iy + 1] = sum[1];
397 ai += incai;
398 iy += incy;
399 } /* end for */
400 } else { /* alpha != 1 */
401 ai = 0;
402 iy = ky;
403 for (i = 0; i < leny; i++) {
404 sum[0] = sum[1] = 0.0;
405 sum2[0] = sum2[1] = 0.0;
406 aij = ai;
407 jx = kx;
408 for (j = 0; j < lenx; j++) {
409 a_elem[0] = a_i[aij];
410 a_elem[1] = a_i[aij + 1];
411
412 x_elem = head_x_i[jx];
413 {
414 prod[0] = a_elem[0] * x_elem;
415 prod[1] = a_elem[1] * x_elem;
416 }
417 sum[0] = sum[0] + prod[0];
418 sum[1] = sum[1] + prod[1];
419 x_elem = tail_x_i[jx];
420 {
421 prod[0] = a_elem[0] * x_elem;
422 prod[1] = a_elem[1] * x_elem;
423 }
424 sum2[0] = sum2[0] + prod[0];
425 sum2[1] = sum2[1] + prod[1];
426 aij += incaij;
427 jx += incx;
428 }
429 {
430 tmp1[0] = sum[0] * alpha_i[0] - sum[1] * alpha_i[1];
431 tmp1[1] = sum[0] * alpha_i[1] + sum[1] * alpha_i[0];
432 }
433
434 {
435 tmp2[0] = sum2[0] * alpha_i[0] - sum2[1] * alpha_i[1];
436 tmp2[1] = sum2[0] * alpha_i[1] + sum2[1] * alpha_i[0];
437 }
438
439 tmp1[0] = tmp1[0] + tmp2[0];
440 tmp1[1] = tmp1[1] + tmp2[1];
441 y_i[iy] = tmp1[0];
442 y_i[iy + 1] = tmp1[1];
443 ai += incai;
444 iy += incy;
445 }
446 }
447 } else { /* beta != 0 */
448 if ((alpha_i[0] == 1.0 && alpha_i[1] == 0.0)) {
449 /* save m multiplies if alpha = 1 */
450 ai = 0;
451 iy = ky;
452 for (i = 0; i < leny; i++) {
453 sum[0] = sum[1] = 0.0;;
454 sum2[0] = sum2[1] = 0.0;;
455 aij = ai;
456 jx = kx;
457 for (j = 0; j < lenx; j++) {
458 a_elem[0] = a_i[aij];
459 a_elem[1] = a_i[aij + 1];
460
461 x_elem = head_x_i[jx];
462 {
463 prod[0] = a_elem[0] * x_elem;
464 prod[1] = a_elem[1] * x_elem;
465 }
466 sum[0] = sum[0] + prod[0];
467 sum[1] = sum[1] + prod[1];
468 x_elem = tail_x_i[jx];
469 {
470 prod[0] = a_elem[0] * x_elem;
471 prod[1] = a_elem[1] * x_elem;
472 }
473 sum2[0] = sum2[0] + prod[0];
474 sum2[1] = sum2[1] + prod[1];
475 aij += incaij;
476 jx += incx;
477 }
478 sum[0] = sum[0] + sum2[0];
479 sum[1] = sum[1] + sum2[1];
480 y_elem[0] = y_i[iy];
481 y_elem[1] = y_i[iy + 1];
482 {
483 tmp1[0] = y_elem[0] * beta_i[0] - y_elem[1] * beta_i[1];
484 tmp1[1] = y_elem[0] * beta_i[1] + y_elem[1] * beta_i[0];
485 }
486
487 tmp2[0] = sum[0] + tmp1[0];
488 tmp2[1] = sum[1] + tmp1[1];
489 y_i[iy] = tmp2[0];
490 y_i[iy + 1] = tmp2[1];
491 ai += incai;
492 iy += incy;
493 }
494 } else { /* alpha != 1, the most general form:
495 y = alpha*A*head_x + alpha*A*tail_x + beta*y */
496 ai = 0;
497 iy = ky;
498 for (i = 0; i < leny; i++) {
499 sum[0] = sum[1] = 0.0;;
500 sum2[0] = sum2[1] = 0.0;;
501 aij = ai;
502 jx = kx;
503 for (j = 0; j < lenx; j++) {
504 a_elem[0] = a_i[aij];
505 a_elem[1] = a_i[aij + 1];
506
507 x_elem = head_x_i[jx];
508 {
509 prod[0] = a_elem[0] * x_elem;
510 prod[1] = a_elem[1] * x_elem;
511 }
512 sum[0] = sum[0] + prod[0];
513 sum[1] = sum[1] + prod[1];
514 x_elem = tail_x_i[jx];
515 {
516 prod[0] = a_elem[0] * x_elem;
517 prod[1] = a_elem[1] * x_elem;
518 }
519 sum2[0] = sum2[0] + prod[0];
520 sum2[1] = sum2[1] + prod[1];
521 aij += incaij;
522 jx += incx;
523 }
524 {
525 tmp1[0] = sum[0] * alpha_i[0] - sum[1] * alpha_i[1];
526 tmp1[1] = sum[0] * alpha_i[1] + sum[1] * alpha_i[0];
527 }
528
529 {
530 tmp2[0] = sum2[0] * alpha_i[0] - sum2[1] * alpha_i[1];
531 tmp2[1] = sum2[0] * alpha_i[1] + sum2[1] * alpha_i[0];
532 }
533
534 tmp1[0] = tmp1[0] + tmp2[0];
535 tmp1[1] = tmp1[1] + tmp2[1];
536 y_elem[0] = y_i[iy];
537 y_elem[1] = y_i[iy + 1];
538 {
539 tmp2[0] = y_elem[0] * beta_i[0] - y_elem[1] * beta_i[1];
540 tmp2[1] = y_elem[0] * beta_i[1] + y_elem[1] * beta_i[0];
541 }
542
543 tmp1[0] = tmp1[0] + tmp2[0];
544 tmp1[1] = tmp1[1] + tmp2[1];
545 y_i[iy] = tmp1[0];
546 y_i[iy + 1] = tmp1[1];
547 ai += incai;
548 iy += incy;
549 }
550 }
551 }
552
553 }
554 }
555
556
557
558 break;
559 }
560 case blas_prec_double:
561 case blas_prec_indigenous:{
562
563 int i, j;
564 int iy, jx, kx, ky;
565 int lenx, leny;
566 int ai, aij;
567 int incai, incaij;
568
569 const float *a_i = (float *) a;
570 const float *head_x_i = head_x;
571 const float *tail_x_i = tail_x;
572 float *y_i = (float *) y;
573 float *alpha_i = (float *) alpha;
574 float *beta_i = (float *) beta;
575 float a_elem[2];
576 float x_elem;
577 float y_elem[2];
578 double prod[2];
579 double sum[2];
580 double sum2[2];
581 double tmp1[2];
582 double tmp2[2];
583
584
585 /* all error calls */
586 if (m < 0)
587 BLAS_error(routine_name, -3, m, 0);
588 else if (n <= 0)
589 BLAS_error(routine_name, -4, n, 0);
590 else if (incx == 0)
591 BLAS_error(routine_name, -10, incx, 0);
592 else if (incy == 0)
593 BLAS_error(routine_name, -13, incy, 0);
594
595 if ((order == blas_rowmajor) && (trans == blas_no_trans)) {
596 lenx = n;
597 leny = m;
598 incai = lda;
599 incaij = 1;
600 } else if ((order == blas_rowmajor) && (trans != blas_no_trans)) {
601 lenx = m;
602 leny = n;
603 incai = 1;
604 incaij = lda;
605 } else if ((order == blas_colmajor) && (trans == blas_no_trans)) {
606 lenx = n;
607 leny = m;
608 incai = 1;
609 incaij = lda;
610 } else { /* colmajor and blas_trans */
611 lenx = m;
612 leny = n;
613 incai = lda;
614 incaij = 1;
615 }
616
617 if (lda < leny)
618 BLAS_error(routine_name, -7, lda, NULL);
619
620
621
622
623 incy *= 2;
624 incai *= 2;
625 incaij *= 2;
626
627 if (incx > 0)
628 kx = 0;
629 else
630 kx = (1 - lenx) * incx;
631 if (incy > 0)
632 ky = 0;
633 else
634 ky = (1 - leny) * incy;
635
636 /* No extra-precision needed for alpha = 0 */
637 if (alpha_i[0] == 0.0 && alpha_i[1] == 0.0) {
638 if (beta_i[0] == 0.0 && beta_i[1] == 0.0) {
639 iy = ky;
640 for (i = 0; i < leny; i++) {
641 y_i[iy] = 0.0;
642 y_i[iy + 1] = 0.0;
643 iy += incy;
644 }
645 } else if (!(beta_i[0] == 0.0 && beta_i[1] == 0.0)) {
646 iy = ky;
647 for (i = 0; i < leny; i++) {
648 y_elem[0] = y_i[iy];
649 y_elem[1] = y_i[iy + 1];
650 {
651 tmp1[0] =
652 (double) y_elem[0] * beta_i[0] -
653 (double) y_elem[1] * beta_i[1];
654 tmp1[1] =
655 (double) y_elem[0] * beta_i[1] +
656 (double) y_elem[1] * beta_i[0];
657 }
658 y_i[iy] = tmp1[0];
659 y_i[iy + 1] = tmp1[1];
660 iy += incy;
661 }
662 }
663 } else { /* alpha != 0 */
664 if (trans == blas_conj_trans) {
665
666 /* if beta = 0, we can save m multiplies:
667 y = alpha*A*head_x + alpha*A*tail_x */
668 if (beta_i[0] == 0.0 && beta_i[1] == 0.0) {
669 if ((alpha_i[0] == 1.0 && alpha_i[1] == 0.0)) {
670 /* save m more multiplies if alpha = 1 */
671 ai = 0;
672 iy = ky;
673 for (i = 0; i < leny; i++) {
674 sum[0] = sum[1] = 0.0;
675 sum2[0] = sum2[1] = 0.0;
676 aij = ai;
677 jx = kx;
678 for (j = 0; j < lenx; j++) {
679 a_elem[0] = a_i[aij];
680 a_elem[1] = a_i[aij + 1];
681 a_elem[1] = -a_elem[1];
682 x_elem = head_x_i[jx];
683 {
684 prod[0] = (double) a_elem[0] * x_elem;
685 prod[1] = (double) a_elem[1] * x_elem;
686 }
687 sum[0] = sum[0] + prod[0];
688 sum[1] = sum[1] + prod[1];
689 x_elem = tail_x_i[jx];
690 {
691 prod[0] = (double) a_elem[0] * x_elem;
692 prod[1] = (double) a_elem[1] * x_elem;
693 }
694 sum2[0] = sum2[0] + prod[0];
695 sum2[1] = sum2[1] + prod[1];
696 aij += incaij;
697 jx += incx;
698 }
699 sum[0] = sum[0] + sum2[0];
700 sum[1] = sum[1] + sum2[1];
701 y_i[iy] = sum[0];
702 y_i[iy + 1] = sum[1];
703 ai += incai;
704 iy += incy;
705 } /* end for */
706 } else { /* alpha != 1 */
707 ai = 0;
708 iy = ky;
709 for (i = 0; i < leny; i++) {
710 sum[0] = sum[1] = 0.0;
711 sum2[0] = sum2[1] = 0.0;
712 aij = ai;
713 jx = kx;
714 for (j = 0; j < lenx; j++) {
715 a_elem[0] = a_i[aij];
716 a_elem[1] = a_i[aij + 1];
717 a_elem[1] = -a_elem[1];
718 x_elem = head_x_i[jx];
719 {
720 prod[0] = (double) a_elem[0] * x_elem;
721 prod[1] = (double) a_elem[1] * x_elem;
722 }
723 sum[0] = sum[0] + prod[0];
724 sum[1] = sum[1] + prod[1];
725 x_elem = tail_x_i[jx];
726 {
727 prod[0] = (double) a_elem[0] * x_elem;
728 prod[1] = (double) a_elem[1] * x_elem;
729 }
730 sum2[0] = sum2[0] + prod[0];
731 sum2[1] = sum2[1] + prod[1];
732 aij += incaij;
733 jx += incx;
734 }
735 {
736 tmp1[0] =
737 (double) sum[0] * alpha_i[0] -
738 (double) sum[1] * alpha_i[1];
739 tmp1[1] =
740 (double) sum[0] * alpha_i[1] +
741 (double) sum[1] * alpha_i[0];
742 }
743 {
744 tmp2[0] =
745 (double) sum2[0] * alpha_i[0] -
746 (double) sum2[1] * alpha_i[1];
747 tmp2[1] =
748 (double) sum2[0] * alpha_i[1] +
749 (double) sum2[1] * alpha_i[0];
750 }
751 tmp1[0] = tmp1[0] + tmp2[0];
752 tmp1[1] = tmp1[1] + tmp2[1];
753 y_i[iy] = tmp1[0];
754 y_i[iy + 1] = tmp1[1];
755 ai += incai;
756 iy += incy;
757 }
758 }
759 } else { /* beta != 0 */
760 if ((alpha_i[0] == 1.0 && alpha_i[1] == 0.0)) {
761 /* save m multiplies if alpha = 1 */
762 ai = 0;
763 iy = ky;
764 for (i = 0; i < leny; i++) {
765 sum[0] = sum[1] = 0.0;;
766 sum2[0] = sum2[1] = 0.0;;
767 aij = ai;
768 jx = kx;
769 for (j = 0; j < lenx; j++) {
770 a_elem[0] = a_i[aij];
771 a_elem[1] = a_i[aij + 1];
772 a_elem[1] = -a_elem[1];
773 x_elem = head_x_i[jx];
774 {
775 prod[0] = (double) a_elem[0] * x_elem;
776 prod[1] = (double) a_elem[1] * x_elem;
777 }
778 sum[0] = sum[0] + prod[0];
779 sum[1] = sum[1] + prod[1];
780 x_elem = tail_x_i[jx];
781 {
782 prod[0] = (double) a_elem[0] * x_elem;
783 prod[1] = (double) a_elem[1] * x_elem;
784 }
785 sum2[0] = sum2[0] + prod[0];
786 sum2[1] = sum2[1] + prod[1];
787 aij += incaij;
788 jx += incx;
789 }
790 sum[0] = sum[0] + sum2[0];
791 sum[1] = sum[1] + sum2[1];
792 y_elem[0] = y_i[iy];
793 y_elem[1] = y_i[iy + 1];
794 {
795 tmp1[0] =
796 (double) y_elem[0] * beta_i[0] -
797 (double) y_elem[1] * beta_i[1];
798 tmp1[1] =
799 (double) y_elem[0] * beta_i[1] +
800 (double) y_elem[1] * beta_i[0];
801 }
802 tmp2[0] = sum[0] + tmp1[0];
803 tmp2[1] = sum[1] + tmp1[1];
804 y_i[iy] = tmp2[0];
805 y_i[iy + 1] = tmp2[1];
806 ai += incai;
807 iy += incy;
808 }
809 } else { /* alpha != 1, the most general form:
810 y = alpha*A*head_x + alpha*A*tail_x + beta*y */
811 ai = 0;
812 iy = ky;
813 for (i = 0; i < leny; i++) {
814 sum[0] = sum[1] = 0.0;;
815 sum2[0] = sum2[1] = 0.0;;
816 aij = ai;
817 jx = kx;
818 for (j = 0; j < lenx; j++) {
819 a_elem[0] = a_i[aij];
820 a_elem[1] = a_i[aij + 1];
821 a_elem[1] = -a_elem[1];
822 x_elem = head_x_i[jx];
823 {
824 prod[0] = (double) a_elem[0] * x_elem;
825 prod[1] = (double) a_elem[1] * x_elem;
826 }
827 sum[0] = sum[0] + prod[0];
828 sum[1] = sum[1] + prod[1];
829 x_elem = tail_x_i[jx];
830 {
831 prod[0] = (double) a_elem[0] * x_elem;
832 prod[1] = (double) a_elem[1] * x_elem;
833 }
834 sum2[0] = sum2[0] + prod[0];
835 sum2[1] = sum2[1] + prod[1];
836 aij += incaij;
837 jx += incx;
838 }
839 {
840 tmp1[0] =
841 (double) sum[0] * alpha_i[0] -
842 (double) sum[1] * alpha_i[1];
843 tmp1[1] =
844 (double) sum[0] * alpha_i[1] +
845 (double) sum[1] * alpha_i[0];
846 }
847 {
848 tmp2[0] =
849 (double) sum2[0] * alpha_i[0] -
850 (double) sum2[1] * alpha_i[1];
851 tmp2[1] =
852 (double) sum2[0] * alpha_i[1] +
853 (double) sum2[1] * alpha_i[0];
854 }
855 tmp1[0] = tmp1[0] + tmp2[0];
856 tmp1[1] = tmp1[1] + tmp2[1];
857 y_elem[0] = y_i[iy];
858 y_elem[1] = y_i[iy + 1];
859 {
860 tmp2[0] =
861 (double) y_elem[0] * beta_i[0] -
862 (double) y_elem[1] * beta_i[1];
863 tmp2[1] =
864 (double) y_elem[0] * beta_i[1] +
865 (double) y_elem[1] * beta_i[0];
866 }
867 tmp1[0] = tmp1[0] + tmp2[0];
868 tmp1[1] = tmp1[1] + tmp2[1];
869 y_i[iy] = tmp1[0];
870 y_i[iy + 1] = tmp1[1];
871 ai += incai;
872 iy += incy;
873 }
874 }
875 }
876
877 } else {
878
879 /* if beta = 0, we can save m multiplies:
880 y = alpha*A*head_x + alpha*A*tail_x */
881 if (beta_i[0] == 0.0 && beta_i[1] == 0.0) {
882 if ((alpha_i[0] == 1.0 && alpha_i[1] == 0.0)) {
883 /* save m more multiplies if alpha = 1 */
884 ai = 0;
885 iy = ky;
886 for (i = 0; i < leny; i++) {
887 sum[0] = sum[1] = 0.0;
888 sum2[0] = sum2[1] = 0.0;
889 aij = ai;
890 jx = kx;
891 for (j = 0; j < lenx; j++) {
892 a_elem[0] = a_i[aij];
893 a_elem[1] = a_i[aij + 1];
894
895 x_elem = head_x_i[jx];
896 {
897 prod[0] = (double) a_elem[0] * x_elem;
898 prod[1] = (double) a_elem[1] * x_elem;
899 }
900 sum[0] = sum[0] + prod[0];
901 sum[1] = sum[1] + prod[1];
902 x_elem = tail_x_i[jx];
903 {
904 prod[0] = (double) a_elem[0] * x_elem;
905 prod[1] = (double) a_elem[1] * x_elem;
906 }
907 sum2[0] = sum2[0] + prod[0];
908 sum2[1] = sum2[1] + prod[1];
909 aij += incaij;
910 jx += incx;
911 }
912 sum[0] = sum[0] + sum2[0];
913 sum[1] = sum[1] + sum2[1];
914 y_i[iy] = sum[0];
915 y_i[iy + 1] = sum[1];
916 ai += incai;
917 iy += incy;
918 } /* end for */
919 } else { /* alpha != 1 */
920 ai = 0;
921 iy = ky;
922 for (i = 0; i < leny; i++) {
923 sum[0] = sum[1] = 0.0;
924 sum2[0] = sum2[1] = 0.0;
925 aij = ai;
926 jx = kx;
927 for (j = 0; j < lenx; j++) {
928 a_elem[0] = a_i[aij];
929 a_elem[1] = a_i[aij + 1];
930
931 x_elem = head_x_i[jx];
932 {
933 prod[0] = (double) a_elem[0] * x_elem;
934 prod[1] = (double) a_elem[1] * x_elem;
935 }
936 sum[0] = sum[0] + prod[0];
937 sum[1] = sum[1] + prod[1];
938 x_elem = tail_x_i[jx];
939 {
940 prod[0] = (double) a_elem[0] * x_elem;
941 prod[1] = (double) a_elem[1] * x_elem;
942 }
943 sum2[0] = sum2[0] + prod[0];
944 sum2[1] = sum2[1] + prod[1];
945 aij += incaij;
946 jx += incx;
947 }
948 {
949 tmp1[0] =
950 (double) sum[0] * alpha_i[0] -
951 (double) sum[1] * alpha_i[1];
952 tmp1[1] =
953 (double) sum[0] * alpha_i[1] +
954 (double) sum[1] * alpha_i[0];
955 }
956 {
957 tmp2[0] =
958 (double) sum2[0] * alpha_i[0] -
959 (double) sum2[1] * alpha_i[1];
960 tmp2[1] =
961 (double) sum2[0] * alpha_i[1] +
962 (double) sum2[1] * alpha_i[0];
963 }
964 tmp1[0] = tmp1[0] + tmp2[0];
965 tmp1[1] = tmp1[1] + tmp2[1];
966 y_i[iy] = tmp1[0];
967 y_i[iy + 1] = tmp1[1];
968 ai += incai;
969 iy += incy;
970 }
971 }
972 } else { /* beta != 0 */
973 if ((alpha_i[0] == 1.0 && alpha_i[1] == 0.0)) {
974 /* save m multiplies if alpha = 1 */
975 ai = 0;
976 iy = ky;
977 for (i = 0; i < leny; i++) {
978 sum[0] = sum[1] = 0.0;;
979 sum2[0] = sum2[1] = 0.0;;
980 aij = ai;
981 jx = kx;
982 for (j = 0; j < lenx; j++) {
983 a_elem[0] = a_i[aij];
984 a_elem[1] = a_i[aij + 1];
985
986 x_elem = head_x_i[jx];
987 {
988 prod[0] = (double) a_elem[0] * x_elem;
989 prod[1] = (double) a_elem[1] * x_elem;
990 }
991 sum[0] = sum[0] + prod[0];
992 sum[1] = sum[1] + prod[1];
993 x_elem = tail_x_i[jx];
994 {
995 prod[0] = (double) a_elem[0] * x_elem;
996 prod[1] = (double) a_elem[1] * x_elem;
997 }
998 sum2[0] = sum2[0] + prod[0];
999 sum2[1] = sum2[1] + prod[1];
1000 aij += incaij;
1001 jx += incx;
1002 }
1003 sum[0] = sum[0] + sum2[0];
1004 sum[1] = sum[1] + sum2[1];
1005 y_elem[0] = y_i[iy];
1006 y_elem[1] = y_i[iy + 1];
1007 {
1008 tmp1[0] =
1009 (double) y_elem[0] * beta_i[0] -
1010 (double) y_elem[1] * beta_i[1];
1011 tmp1[1] =
1012 (double) y_elem[0] * beta_i[1] +
1013 (double) y_elem[1] * beta_i[0];
1014 }
1015 tmp2[0] = sum[0] + tmp1[0];
1016 tmp2[1] = sum[1] + tmp1[1];
1017 y_i[iy] = tmp2[0];
1018 y_i[iy + 1] = tmp2[1];
1019 ai += incai;
1020 iy += incy;
1021 }
1022 } else { /* alpha != 1, the most general form:
1023 y = alpha*A*head_x + alpha*A*tail_x + beta*y */
1024 ai = 0;
1025 iy = ky;
1026 for (i = 0; i < leny; i++) {
1027 sum[0] = sum[1] = 0.0;;
1028 sum2[0] = sum2[1] = 0.0;;
1029 aij = ai;
1030 jx = kx;
1031 for (j = 0; j < lenx; j++) {
1032 a_elem[0] = a_i[aij];
1033 a_elem[1] = a_i[aij + 1];
1034
1035 x_elem = head_x_i[jx];
1036 {
1037 prod[0] = (double) a_elem[0] * x_elem;
1038 prod[1] = (double) a_elem[1] * x_elem;
1039 }
1040 sum[0] = sum[0] + prod[0];
1041 sum[1] = sum[1] + prod[1];
1042 x_elem = tail_x_i[jx];
1043 {
1044 prod[0] = (double) a_elem[0] * x_elem;
1045 prod[1] = (double) a_elem[1] * x_elem;
1046 }
1047 sum2[0] = sum2[0] + prod[0];
1048 sum2[1] = sum2[1] + prod[1];
1049 aij += incaij;
1050 jx += incx;
1051 }
1052 {
1053 tmp1[0] =
1054 (double) sum[0] * alpha_i[0] -
1055 (double) sum[1] * alpha_i[1];
1056 tmp1[1] =
1057 (double) sum[0] * alpha_i[1] +
1058 (double) sum[1] * alpha_i[0];
1059 }
1060 {
1061 tmp2[0] =
1062 (double) sum2[0] * alpha_i[0] -
1063 (double) sum2[1] * alpha_i[1];
1064 tmp2[1] =
1065 (double) sum2[0] * alpha_i[1] +
1066 (double) sum2[1] * alpha_i[0];
1067 }
1068 tmp1[0] = tmp1[0] + tmp2[0];
1069 tmp1[1] = tmp1[1] + tmp2[1];
1070 y_elem[0] = y_i[iy];
1071 y_elem[1] = y_i[iy + 1];
1072 {
1073 tmp2[0] =
1074 (double) y_elem[0] * beta_i[0] -
1075 (double) y_elem[1] * beta_i[1];
1076 tmp2[1] =
1077 (double) y_elem[0] * beta_i[1] +
1078 (double) y_elem[1] * beta_i[0];
1079 }
1080 tmp1[0] = tmp1[0] + tmp2[0];
1081 tmp1[1] = tmp1[1] + tmp2[1];
1082 y_i[iy] = tmp1[0];
1083 y_i[iy + 1] = tmp1[1];
1084 ai += incai;
1085 iy += incy;
1086 }
1087 }
1088 }
1089
1090 }
1091 }
1092
1093
1094
1095 break;
1096 }
1097 case blas_prec_extra:{
1098
1099 int i, j;
1100 int iy, jx, kx, ky;
1101 int lenx, leny;
1102 int ai, aij;
1103 int incai, incaij;
1104
1105 const float *a_i = (float *) a;
1106 const float *head_x_i = head_x;
1107 const float *tail_x_i = tail_x;
1108 float *y_i = (float *) y;
1109 float *alpha_i = (float *) alpha;
1110 float *beta_i = (float *) beta;
1111 float a_elem[2];
1112 float x_elem;
1113 float y_elem[2];
1114 double head_prod[2], tail_prod[2];
1115 double head_sum[2], tail_sum[2];
1116 double head_sum2[2], tail_sum2[2];
1117 double head_tmp1[2], tail_tmp1[2];
1118 double head_tmp2[2], tail_tmp2[2];
1119 FPU_FIX_DECL;
1120
1121 /* all error calls */
1122 if (m < 0)
1123 BLAS_error(routine_name, -3, m, 0);
1124 else if (n <= 0)
1125 BLAS_error(routine_name, -4, n, 0);
1126 else if (incx == 0)
1127 BLAS_error(routine_name, -10, incx, 0);
1128 else if (incy == 0)
1129 BLAS_error(routine_name, -13, incy, 0);
1130
1131 if ((order == blas_rowmajor) && (trans == blas_no_trans)) {
1132 lenx = n;
1133 leny = m;
1134 incai = lda;
1135 incaij = 1;
1136 } else if ((order == blas_rowmajor) && (trans != blas_no_trans)) {
1137 lenx = m;
1138 leny = n;
1139 incai = 1;
1140 incaij = lda;
1141 } else if ((order == blas_colmajor) && (trans == blas_no_trans)) {
1142 lenx = n;
1143 leny = m;
1144 incai = 1;
1145 incaij = lda;
1146 } else { /* colmajor and blas_trans */
1147 lenx = m;
1148 leny = n;
1149 incai = lda;
1150 incaij = 1;
1151 }
1152
1153 if (lda < leny)
1154 BLAS_error(routine_name, -7, lda, NULL);
1155
1156 FPU_FIX_START;
1157
1158
1159 incy *= 2;
1160 incai *= 2;
1161 incaij *= 2;
1162
1163 if (incx > 0)
1164 kx = 0;
1165 else
1166 kx = (1 - lenx) * incx;
1167 if (incy > 0)
1168 ky = 0;
1169 else
1170 ky = (1 - leny) * incy;
1171
1172 /* No extra-precision needed for alpha = 0 */
1173 if (alpha_i[0] == 0.0 && alpha_i[1] == 0.0) {
1174 if (beta_i[0] == 0.0 && beta_i[1] == 0.0) {
1175 iy = ky;
1176 for (i = 0; i < leny; i++) {
1177 y_i[iy] = 0.0;
1178 y_i[iy + 1] = 0.0;
1179 iy += incy;
1180 }
1181 } else if (!(beta_i[0] == 0.0 && beta_i[1] == 0.0)) {
1182 iy = ky;
1183 for (i = 0; i < leny; i++) {
1184 y_elem[0] = y_i[iy];
1185 y_elem[1] = y_i[iy + 1];
1186 {
1187 double head_e1, tail_e1;
1188 double d1;
1189 double d2;
1190 /* Real part */
1191 d1 = (double) y_elem[0] * beta_i[0];
1192 d2 = (double) -y_elem[1] * beta_i[1];
1193 {
1194 /* Compute double-double = double + double. */
1195 double e, t1, t2;
1196
1197 /* Knuth trick. */
1198 t1 = d1 + d2;
1199 e = t1 - d1;
1200 t2 = ((d2 - e) + (d1 - (t1 - e)));
1201
1202 /* The result is t1 + t2, after normalization. */
1203 head_e1 = t1 + t2;
1204 tail_e1 = t2 - (head_e1 - t1);
1205 }
1206 head_tmp1[0] = head_e1;
1207 tail_tmp1[0] = tail_e1;
1208 /* imaginary part */
1209 d1 = (double) y_elem[0] * beta_i[1];
1210 d2 = (double) y_elem[1] * beta_i[0];
1211 {
1212 /* Compute double-double = double + double. */
1213 double e, t1, t2;
1214
1215 /* Knuth trick. */
1216 t1 = d1 + d2;
1217 e = t1 - d1;
1218 t2 = ((d2 - e) + (d1 - (t1 - e)));
1219
1220 /* The result is t1 + t2, after normalization. */
1221 head_e1 = t1 + t2;
1222 tail_e1 = t2 - (head_e1 - t1);
1223 }
1224 head_tmp1[1] = head_e1;
1225 tail_tmp1[1] = tail_e1;
1226 }
1227 y_i[iy] = head_tmp1[0];
1228 y_i[iy + 1] = head_tmp1[1];
1229 iy += incy;
1230 }
1231 }
1232 } else { /* alpha != 0 */
1233 if (trans == blas_conj_trans) {
1234
1235 /* if beta = 0, we can save m multiplies:
1236 y = alpha*A*head_x + alpha*A*tail_x */
1237 if (beta_i[0] == 0.0 && beta_i[1] == 0.0) {
1238 if ((alpha_i[0] == 1.0 && alpha_i[1] == 0.0)) {
1239 /* save m more multiplies if alpha = 1 */
1240 ai = 0;
1241 iy = ky;
1242 for (i = 0; i < leny; i++) {
1243 head_sum[0] = head_sum[1] = tail_sum[0] = tail_sum[1] = 0.0;
1244 head_sum2[0] = head_sum2[1] = tail_sum2[0] = tail_sum2[1] =
1245 0.0;
1246 aij = ai;
1247 jx = kx;
1248 for (j = 0; j < lenx; j++) {
1249 a_elem[0] = a_i[aij];
1250 a_elem[1] = a_i[aij + 1];
1251 a_elem[1] = -a_elem[1];
1252 x_elem = head_x_i[jx];
1253 {
1254 head_prod[0] = (double) a_elem[0] * x_elem;
1255 tail_prod[0] = 0.0;
1256 head_prod[1] = (double) a_elem[1] * x_elem;
1257 tail_prod[1] = 0.0;
1258 }
1259 {
1260 double head_t, tail_t;
1261 double head_a, tail_a;
1262 double head_b, tail_b;
1263 /* Real part */
1264 head_a = head_sum[0];
1265 tail_a = tail_sum[0];
1266 head_b = head_prod[0];
1267 tail_b = tail_prod[0];
1268 {
1269 /* Compute double-double = double-double + double-double. */
1270 double bv;
1271 double s1, s2, t1, t2;
1272
1273 /* Add two hi words. */
1274 s1 = head_a + head_b;
1275 bv = s1 - head_a;
1276 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1277
1278 /* Add two lo words. */
1279 t1 = tail_a + tail_b;
1280 bv = t1 - tail_a;
1281 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1282
1283 s2 += t1;
1284
1285 /* Renormalize (s1, s2) to (t1, s2) */
1286 t1 = s1 + s2;
1287 s2 = s2 - (t1 - s1);
1288
1289 t2 += s2;
1290
1291 /* Renormalize (t1, t2) */
1292 head_t = t1 + t2;
1293 tail_t = t2 - (head_t - t1);
1294 }
1295 head_sum[0] = head_t;
1296 tail_sum[0] = tail_t;
1297 /* Imaginary part */
1298 head_a = head_sum[1];
1299 tail_a = tail_sum[1];
1300 head_b = head_prod[1];
1301 tail_b = tail_prod[1];
1302 {
1303 /* Compute double-double = double-double + double-double. */
1304 double bv;
1305 double s1, s2, t1, t2;
1306
1307 /* Add two hi words. */
1308 s1 = head_a + head_b;
1309 bv = s1 - head_a;
1310 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1311
1312 /* Add two lo words. */
1313 t1 = tail_a + tail_b;
1314 bv = t1 - tail_a;
1315 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1316
1317 s2 += t1;
1318
1319 /* Renormalize (s1, s2) to (t1, s2) */
1320 t1 = s1 + s2;
1321 s2 = s2 - (t1 - s1);
1322
1323 t2 += s2;
1324
1325 /* Renormalize (t1, t2) */
1326 head_t = t1 + t2;
1327 tail_t = t2 - (head_t - t1);
1328 }
1329 head_sum[1] = head_t;
1330 tail_sum[1] = tail_t;
1331 }
1332 x_elem = tail_x_i[jx];
1333 {
1334 head_prod[0] = (double) a_elem[0] * x_elem;
1335 tail_prod[0] = 0.0;
1336 head_prod[1] = (double) a_elem[1] * x_elem;
1337 tail_prod[1] = 0.0;
1338 }
1339 {
1340 double head_t, tail_t;
1341 double head_a, tail_a;
1342 double head_b, tail_b;
1343 /* Real part */
1344 head_a = head_sum2[0];
1345 tail_a = tail_sum2[0];
1346 head_b = head_prod[0];
1347 tail_b = tail_prod[0];
1348 {
1349 /* Compute double-double = double-double + double-double. */
1350 double bv;
1351 double s1, s2, t1, t2;
1352
1353 /* Add two hi words. */
1354 s1 = head_a + head_b;
1355 bv = s1 - head_a;
1356 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1357
1358 /* Add two lo words. */
1359 t1 = tail_a + tail_b;
1360 bv = t1 - tail_a;
1361 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1362
1363 s2 += t1;
1364
1365 /* Renormalize (s1, s2) to (t1, s2) */
1366 t1 = s1 + s2;
1367 s2 = s2 - (t1 - s1);
1368
1369 t2 += s2;
1370
1371 /* Renormalize (t1, t2) */
1372 head_t = t1 + t2;
1373 tail_t = t2 - (head_t - t1);
1374 }
1375 head_sum2[0] = head_t;
1376 tail_sum2[0] = tail_t;
1377 /* Imaginary part */
1378 head_a = head_sum2[1];
1379 tail_a = tail_sum2[1];
1380 head_b = head_prod[1];
1381 tail_b = tail_prod[1];
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_a + head_b;
1389 bv = s1 - head_a;
1390 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1391
1392 /* Add two lo words. */
1393 t1 = tail_a + tail_b;
1394 bv = t1 - tail_a;
1395 t2 = ((tail_b - bv) + (tail_a - (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_t = t1 + t2;
1407 tail_t = t2 - (head_t - t1);
1408 }
1409 head_sum2[1] = head_t;
1410 tail_sum2[1] = tail_t;
1411 }
1412 aij += incaij;
1413 jx += incx;
1414 }
1415 {
1416 double head_t, tail_t;
1417 double head_a, tail_a;
1418 double head_b, tail_b;
1419 /* Real part */
1420 head_a = head_sum[0];
1421 tail_a = tail_sum[0];
1422 head_b = head_sum2[0];
1423 tail_b = tail_sum2[0];
1424 {
1425 /* Compute double-double = double-double + double-double. */
1426 double bv;
1427 double s1, s2, t1, t2;
1428
1429 /* Add two hi words. */
1430 s1 = head_a + head_b;
1431 bv = s1 - head_a;
1432 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1433
1434 /* Add two lo words. */
1435 t1 = tail_a + tail_b;
1436 bv = t1 - tail_a;
1437 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1438
1439 s2 += t1;
1440
1441 /* Renormalize (s1, s2) to (t1, s2) */
1442 t1 = s1 + s2;
1443 s2 = s2 - (t1 - s1);
1444
1445 t2 += s2;
1446
1447 /* Renormalize (t1, t2) */
1448 head_t = t1 + t2;
1449 tail_t = t2 - (head_t - t1);
1450 }
1451 head_sum[0] = head_t;
1452 tail_sum[0] = tail_t;
1453 /* Imaginary part */
1454 head_a = head_sum[1];
1455 tail_a = tail_sum[1];
1456 head_b = head_sum2[1];
1457 tail_b = tail_sum2[1];
1458 {
1459 /* Compute double-double = double-double + double-double. */
1460 double bv;
1461 double s1, s2, t1, t2;
1462
1463 /* Add two hi words. */
1464 s1 = head_a + head_b;
1465 bv = s1 - head_a;
1466 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1467
1468 /* Add two lo words. */
1469 t1 = tail_a + tail_b;
1470 bv = t1 - tail_a;
1471 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1472
1473 s2 += t1;
1474
1475 /* Renormalize (s1, s2) to (t1, s2) */
1476 t1 = s1 + s2;
1477 s2 = s2 - (t1 - s1);
1478
1479 t2 += s2;
1480
1481 /* Renormalize (t1, t2) */
1482 head_t = t1 + t2;
1483 tail_t = t2 - (head_t - t1);
1484 }
1485 head_sum[1] = head_t;
1486 tail_sum[1] = tail_t;
1487 }
1488 y_i[iy] = head_sum[0];
1489 y_i[iy + 1] = head_sum[1];
1490 ai += incai;
1491 iy += incy;
1492 } /* end for */
1493 } else { /* alpha != 1 */
1494 ai = 0;
1495 iy = ky;
1496 for (i = 0; i < leny; i++) {
1497 head_sum[0] = head_sum[1] = tail_sum[0] = tail_sum[1] = 0.0;
1498 head_sum2[0] = head_sum2[1] = tail_sum2[0] = tail_sum2[1] =
1499 0.0;
1500 aij = ai;
1501 jx = kx;
1502 for (j = 0; j < lenx; j++) {
1503 a_elem[0] = a_i[aij];
1504 a_elem[1] = a_i[aij + 1];
1505 a_elem[1] = -a_elem[1];
1506 x_elem = head_x_i[jx];
1507 {
1508 head_prod[0] = (double) a_elem[0] * x_elem;
1509 tail_prod[0] = 0.0;
1510 head_prod[1] = (double) a_elem[1] * x_elem;
1511 tail_prod[1] = 0.0;
1512 }
1513 {
1514 double head_t, tail_t;
1515 double head_a, tail_a;
1516 double head_b, tail_b;
1517 /* Real part */
1518 head_a = head_sum[0];
1519 tail_a = tail_sum[0];
1520 head_b = head_prod[0];
1521 tail_b = tail_prod[0];
1522 {
1523 /* Compute double-double = double-double + double-double. */
1524 double bv;
1525 double s1, s2, t1, t2;
1526
1527 /* Add two hi words. */
1528 s1 = head_a + head_b;
1529 bv = s1 - head_a;
1530 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1531
1532 /* Add two lo words. */
1533 t1 = tail_a + tail_b;
1534 bv = t1 - tail_a;
1535 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1536
1537 s2 += t1;
1538
1539 /* Renormalize (s1, s2) to (t1, s2) */
1540 t1 = s1 + s2;
1541 s2 = s2 - (t1 - s1);
1542
1543 t2 += s2;
1544
1545 /* Renormalize (t1, t2) */
1546 head_t = t1 + t2;
1547 tail_t = t2 - (head_t - t1);
1548 }
1549 head_sum[0] = head_t;
1550 tail_sum[0] = tail_t;
1551 /* Imaginary part */
1552 head_a = head_sum[1];
1553 tail_a = tail_sum[1];
1554 head_b = head_prod[1];
1555 tail_b = tail_prod[1];
1556 {
1557 /* Compute double-double = double-double + double-double. */
1558 double bv;
1559 double s1, s2, t1, t2;
1560
1561 /* Add two hi words. */
1562 s1 = head_a + head_b;
1563 bv = s1 - head_a;
1564 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1565
1566 /* Add two lo words. */
1567 t1 = tail_a + tail_b;
1568 bv = t1 - tail_a;
1569 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1570
1571 s2 += t1;
1572
1573 /* Renormalize (s1, s2) to (t1, s2) */
1574 t1 = s1 + s2;
1575 s2 = s2 - (t1 - s1);
1576
1577 t2 += s2;
1578
1579 /* Renormalize (t1, t2) */
1580 head_t = t1 + t2;
1581 tail_t = t2 - (head_t - t1);
1582 }
1583 head_sum[1] = head_t;
1584 tail_sum[1] = tail_t;
1585 }
1586 x_elem = tail_x_i[jx];
1587 {
1588 head_prod[0] = (double) a_elem[0] * x_elem;
1589 tail_prod[0] = 0.0;
1590 head_prod[1] = (double) a_elem[1] * x_elem;
1591 tail_prod[1] = 0.0;
1592 }
1593 {
1594 double head_t, tail_t;
1595 double head_a, tail_a;
1596 double head_b, tail_b;
1597 /* Real part */
1598 head_a = head_sum2[0];
1599 tail_a = tail_sum2[0];
1600 head_b = head_prod[0];
1601 tail_b = tail_prod[0];
1602 {
1603 /* Compute double-double = double-double + double-double. */
1604 double bv;
1605 double s1, s2, t1, t2;
1606
1607 /* Add two hi words. */
1608 s1 = head_a + head_b;
1609 bv = s1 - head_a;
1610 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1611
1612 /* Add two lo words. */
1613 t1 = tail_a + tail_b;
1614 bv = t1 - tail_a;
1615 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1616
1617 s2 += t1;
1618
1619 /* Renormalize (s1, s2) to (t1, s2) */
1620 t1 = s1 + s2;
1621 s2 = s2 - (t1 - s1);
1622
1623 t2 += s2;
1624
1625 /* Renormalize (t1, t2) */
1626 head_t = t1 + t2;
1627 tail_t = t2 - (head_t - t1);
1628 }
1629 head_sum2[0] = head_t;
1630 tail_sum2[0] = tail_t;
1631 /* Imaginary part */
1632 head_a = head_sum2[1];
1633 tail_a = tail_sum2[1];
1634 head_b = head_prod[1];
1635 tail_b = tail_prod[1];
1636 {
1637 /* Compute double-double = double-double + double-double. */
1638 double bv;
1639 double s1, s2, t1, t2;
1640
1641 /* Add two hi words. */
1642 s1 = head_a + head_b;
1643 bv = s1 - head_a;
1644 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1645
1646 /* Add two lo words. */
1647 t1 = tail_a + tail_b;
1648 bv = t1 - tail_a;
1649 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1650
1651 s2 += t1;
1652
1653 /* Renormalize (s1, s2) to (t1, s2) */
1654 t1 = s1 + s2;
1655 s2 = s2 - (t1 - s1);
1656
1657 t2 += s2;
1658
1659 /* Renormalize (t1, t2) */
1660 head_t = t1 + t2;
1661 tail_t = t2 - (head_t - t1);
1662 }
1663 head_sum2[1] = head_t;
1664 tail_sum2[1] = tail_t;
1665 }
1666 aij += incaij;
1667 jx += incx;
1668 }
1669 {
1670 double cd[2];
1671 cd[0] = (double) alpha_i[0];
1672 cd[1] = (double) alpha_i[1];
1673 {
1674 /* Compute complex-extra = complex-extra * complex-double. */
1675 double head_a0, tail_a0;
1676 double head_a1, tail_a1;
1677 double head_t1, tail_t1;
1678 double head_t2, tail_t2;
1679 head_a0 = head_sum[0];
1680 tail_a0 = tail_sum[0];
1681 head_a1 = head_sum[1];
1682 tail_a1 = tail_sum[1];
1683 /* real part */
1684 {
1685 /* Compute double-double = double-double * double. */
1686 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
1687
1688 con = head_a0 * split;
1689 a11 = con - head_a0;
1690 a11 = con - a11;
1691 a21 = head_a0 - a11;
1692 con = cd[0] * split;
1693 b1 = con - cd[0];
1694 b1 = con - b1;
1695 b2 = cd[0] - b1;
1696
1697 c11 = head_a0 * cd[0];
1698 c21 =
1699 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
1700
1701 c2 = tail_a0 * cd[0];
1702 t1 = c11 + c2;
1703 t2 = (c2 - (t1 - c11)) + c21;
1704
1705 head_t1 = t1 + t2;
1706 tail_t1 = t2 - (head_t1 - t1);
1707 }
1708 {
1709 /* Compute double-double = double-double * double. */
1710 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
1711
1712 con = head_a1 * split;
1713 a11 = con - head_a1;
1714 a11 = con - a11;
1715 a21 = head_a1 - a11;
1716 con = cd[1] * split;
1717 b1 = con - cd[1];
1718 b1 = con - b1;
1719 b2 = cd[1] - b1;
1720
1721 c11 = head_a1 * cd[1];
1722 c21 =
1723 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
1724
1725 c2 = tail_a1 * cd[1];
1726 t1 = c11 + c2;
1727 t2 = (c2 - (t1 - c11)) + c21;
1728
1729 head_t2 = t1 + t2;
1730 tail_t2 = t2 - (head_t2 - t1);
1731 }
1732 head_t2 = -head_t2;
1733 tail_t2 = -tail_t2;
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_t1 + head_t2;
1741 bv = s1 - head_t1;
1742 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
1743
1744 /* Add two lo words. */
1745 t1 = tail_t1 + tail_t2;
1746 bv = t1 - tail_t1;
1747 t2 = ((tail_t2 - bv) + (tail_t1 - (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_t1 = t1 + t2;
1759 tail_t1 = t2 - (head_t1 - t1);
1760 }
1761 head_tmp1[0] = head_t1;
1762 tail_tmp1[0] = tail_t1;
1763 /* imaginary part */
1764 {
1765 /* Compute double-double = double-double * double. */
1766 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
1767
1768 con = head_a1 * split;
1769 a11 = con - head_a1;
1770 a11 = con - a11;
1771 a21 = head_a1 - a11;
1772 con = cd[0] * split;
1773 b1 = con - cd[0];
1774 b1 = con - b1;
1775 b2 = cd[0] - b1;
1776
1777 c11 = head_a1 * cd[0];
1778 c21 =
1779 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
1780
1781 c2 = tail_a1 * cd[0];
1782 t1 = c11 + c2;
1783 t2 = (c2 - (t1 - c11)) + c21;
1784
1785 head_t1 = t1 + t2;
1786 tail_t1 = t2 - (head_t1 - t1);
1787 }
1788 {
1789 /* Compute double-double = double-double * double. */
1790 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
1791
1792 con = head_a0 * split;
1793 a11 = con - head_a0;
1794 a11 = con - a11;
1795 a21 = head_a0 - a11;
1796 con = cd[1] * split;
1797 b1 = con - cd[1];
1798 b1 = con - b1;
1799 b2 = cd[1] - b1;
1800
1801 c11 = head_a0 * cd[1];
1802 c21 =
1803 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
1804
1805 c2 = tail_a0 * cd[1];
1806 t1 = c11 + c2;
1807 t2 = (c2 - (t1 - c11)) + c21;
1808
1809 head_t2 = t1 + t2;
1810 tail_t2 = t2 - (head_t2 - t1);
1811 }
1812 {
1813 /* Compute double-double = double-double + double-double. */
1814 double bv;
1815 double s1, s2, t1, t2;
1816
1817 /* Add two hi words. */
1818 s1 = head_t1 + head_t2;
1819 bv = s1 - head_t1;
1820 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
1821
1822 /* Add two lo words. */
1823 t1 = tail_t1 + tail_t2;
1824 bv = t1 - tail_t1;
1825 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
1826
1827 s2 += t1;
1828
1829 /* Renormalize (s1, s2) to (t1, s2) */
1830 t1 = s1 + s2;
1831 s2 = s2 - (t1 - s1);
1832
1833 t2 += s2;
1834
1835 /* Renormalize (t1, t2) */
1836 head_t1 = t1 + t2;
1837 tail_t1 = t2 - (head_t1 - t1);
1838 }
1839 head_tmp1[1] = head_t1;
1840 tail_tmp1[1] = tail_t1;
1841 }
1842
1843 }
1844 {
1845 double cd[2];
1846 cd[0] = (double) alpha_i[0];
1847 cd[1] = (double) alpha_i[1];
1848 {
1849 /* Compute complex-extra = complex-extra * complex-double. */
1850 double head_a0, tail_a0;
1851 double head_a1, tail_a1;
1852 double head_t1, tail_t1;
1853 double head_t2, tail_t2;
1854 head_a0 = head_sum2[0];
1855 tail_a0 = tail_sum2[0];
1856 head_a1 = head_sum2[1];
1857 tail_a1 = tail_sum2[1];
1858 /* real part */
1859 {
1860 /* Compute double-double = double-double * double. */
1861 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
1862
1863 con = head_a0 * split;
1864 a11 = con - head_a0;
1865 a11 = con - a11;
1866 a21 = head_a0 - a11;
1867 con = cd[0] * split;
1868 b1 = con - cd[0];
1869 b1 = con - b1;
1870 b2 = cd[0] - b1;
1871
1872 c11 = head_a0 * cd[0];
1873 c21 =
1874 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
1875
1876 c2 = tail_a0 * cd[0];
1877 t1 = c11 + c2;
1878 t2 = (c2 - (t1 - c11)) + c21;
1879
1880 head_t1 = t1 + t2;
1881 tail_t1 = t2 - (head_t1 - t1);
1882 }
1883 {
1884 /* Compute double-double = double-double * double. */
1885 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
1886
1887 con = head_a1 * split;
1888 a11 = con - head_a1;
1889 a11 = con - a11;
1890 a21 = head_a1 - a11;
1891 con = cd[1] * split;
1892 b1 = con - cd[1];
1893 b1 = con - b1;
1894 b2 = cd[1] - b1;
1895
1896 c11 = head_a1 * cd[1];
1897 c21 =
1898 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
1899
1900 c2 = tail_a1 * cd[1];
1901 t1 = c11 + c2;
1902 t2 = (c2 - (t1 - c11)) + c21;
1903
1904 head_t2 = t1 + t2;
1905 tail_t2 = t2 - (head_t2 - t1);
1906 }
1907 head_t2 = -head_t2;
1908 tail_t2 = -tail_t2;
1909 {
1910 /* Compute double-double = double-double + double-double. */
1911 double bv;
1912 double s1, s2, t1, t2;
1913
1914 /* Add two hi words. */
1915 s1 = head_t1 + head_t2;
1916 bv = s1 - head_t1;
1917 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
1918
1919 /* Add two lo words. */
1920 t1 = tail_t1 + tail_t2;
1921 bv = t1 - tail_t1;
1922 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
1923
1924 s2 += t1;
1925
1926 /* Renormalize (s1, s2) to (t1, s2) */
1927 t1 = s1 + s2;
1928 s2 = s2 - (t1 - s1);
1929
1930 t2 += s2;
1931
1932 /* Renormalize (t1, t2) */
1933 head_t1 = t1 + t2;
1934 tail_t1 = t2 - (head_t1 - t1);
1935 }
1936 head_tmp2[0] = head_t1;
1937 tail_tmp2[0] = tail_t1;
1938 /* imaginary part */
1939 {
1940 /* Compute double-double = double-double * double. */
1941 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
1942
1943 con = head_a1 * split;
1944 a11 = con - head_a1;
1945 a11 = con - a11;
1946 a21 = head_a1 - a11;
1947 con = cd[0] * split;
1948 b1 = con - cd[0];
1949 b1 = con - b1;
1950 b2 = cd[0] - b1;
1951
1952 c11 = head_a1 * cd[0];
1953 c21 =
1954 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
1955
1956 c2 = tail_a1 * cd[0];
1957 t1 = c11 + c2;
1958 t2 = (c2 - (t1 - c11)) + c21;
1959
1960 head_t1 = t1 + t2;
1961 tail_t1 = t2 - (head_t1 - t1);
1962 }
1963 {
1964 /* Compute double-double = double-double * double. */
1965 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
1966
1967 con = head_a0 * split;
1968 a11 = con - head_a0;
1969 a11 = con - a11;
1970 a21 = head_a0 - a11;
1971 con = cd[1] * split;
1972 b1 = con - cd[1];
1973 b1 = con - b1;
1974 b2 = cd[1] - b1;
1975
1976 c11 = head_a0 * cd[1];
1977 c21 =
1978 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
1979
1980 c2 = tail_a0 * cd[1];
1981 t1 = c11 + c2;
1982 t2 = (c2 - (t1 - c11)) + c21;
1983
1984 head_t2 = t1 + t2;
1985 tail_t2 = t2 - (head_t2 - t1);
1986 }
1987 {
1988 /* Compute double-double = double-double + double-double. */
1989 double bv;
1990 double s1, s2, t1, t2;
1991
1992 /* Add two hi words. */
1993 s1 = head_t1 + head_t2;
1994 bv = s1 - head_t1;
1995 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
1996
1997 /* Add two lo words. */
1998 t1 = tail_t1 + tail_t2;
1999 bv = t1 - tail_t1;
2000 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
2001
2002 s2 += t1;
2003
2004 /* Renormalize (s1, s2) to (t1, s2) */
2005 t1 = s1 + s2;
2006 s2 = s2 - (t1 - s1);
2007
2008 t2 += s2;
2009
2010 /* Renormalize (t1, t2) */
2011 head_t1 = t1 + t2;
2012 tail_t1 = t2 - (head_t1 - t1);
2013 }
2014 head_tmp2[1] = head_t1;
2015 tail_tmp2[1] = tail_t1;
2016 }
2017
2018 }
2019 {
2020 double head_t, tail_t;
2021 double head_a, tail_a;
2022 double head_b, tail_b;
2023 /* Real part */
2024 head_a = head_tmp1[0];
2025 tail_a = tail_tmp1[0];
2026 head_b = head_tmp2[0];
2027 tail_b = tail_tmp2[0];
2028 {
2029 /* Compute double-double = double-double + double-double. */
2030 double bv;
2031 double s1, s2, t1, t2;
2032
2033 /* Add two hi words. */
2034 s1 = head_a + head_b;
2035 bv = s1 - head_a;
2036 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2037
2038 /* Add two lo words. */
2039 t1 = tail_a + tail_b;
2040 bv = t1 - tail_a;
2041 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2042
2043 s2 += t1;
2044
2045 /* Renormalize (s1, s2) to (t1, s2) */
2046 t1 = s1 + s2;
2047 s2 = s2 - (t1 - s1);
2048
2049 t2 += s2;
2050
2051 /* Renormalize (t1, t2) */
2052 head_t = t1 + t2;
2053 tail_t = t2 - (head_t - t1);
2054 }
2055 head_tmp1[0] = head_t;
2056 tail_tmp1[0] = tail_t;
2057 /* Imaginary part */
2058 head_a = head_tmp1[1];
2059 tail_a = tail_tmp1[1];
2060 head_b = head_tmp2[1];
2061 tail_b = tail_tmp2[1];
2062 {
2063 /* Compute double-double = double-double + double-double. */
2064 double bv;
2065 double s1, s2, t1, t2;
2066
2067 /* Add two hi words. */
2068 s1 = head_a + head_b;
2069 bv = s1 - head_a;
2070 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2071
2072 /* Add two lo words. */
2073 t1 = tail_a + tail_b;
2074 bv = t1 - tail_a;
2075 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2076
2077 s2 += t1;
2078
2079 /* Renormalize (s1, s2) to (t1, s2) */
2080 t1 = s1 + s2;
2081 s2 = s2 - (t1 - s1);
2082
2083 t2 += s2;
2084
2085 /* Renormalize (t1, t2) */
2086 head_t = t1 + t2;
2087 tail_t = t2 - (head_t - t1);
2088 }
2089 head_tmp1[1] = head_t;
2090 tail_tmp1[1] = tail_t;
2091 }
2092 y_i[iy] = head_tmp1[0];
2093 y_i[iy + 1] = head_tmp1[1];
2094 ai += incai;
2095 iy += incy;
2096 }
2097 }
2098 } else { /* beta != 0 */
2099 if ((alpha_i[0] == 1.0 && alpha_i[1] == 0.0)) {
2100 /* save m multiplies if alpha = 1 */
2101 ai = 0;
2102 iy = ky;
2103 for (i = 0; i < leny; i++) {
2104 head_sum[0] = head_sum[1] = tail_sum[0] = tail_sum[1] = 0.0;;
2105 head_sum2[0] = head_sum2[1] = tail_sum2[0] = tail_sum2[1] =
2106 0.0;;
2107 aij = ai;
2108 jx = kx;
2109 for (j = 0; j < lenx; j++) {
2110 a_elem[0] = a_i[aij];
2111 a_elem[1] = a_i[aij + 1];
2112 a_elem[1] = -a_elem[1];
2113 x_elem = head_x_i[jx];
2114 {
2115 head_prod[0] = (double) a_elem[0] * x_elem;
2116 tail_prod[0] = 0.0;
2117 head_prod[1] = (double) a_elem[1] * x_elem;
2118 tail_prod[1] = 0.0;
2119 }
2120 {
2121 double head_t, tail_t;
2122 double head_a, tail_a;
2123 double head_b, tail_b;
2124 /* Real part */
2125 head_a = head_sum[0];
2126 tail_a = tail_sum[0];
2127 head_b = head_prod[0];
2128 tail_b = tail_prod[0];
2129 {
2130 /* Compute double-double = double-double + double-double. */
2131 double bv;
2132 double s1, s2, t1, t2;
2133
2134 /* Add two hi words. */
2135 s1 = head_a + head_b;
2136 bv = s1 - head_a;
2137 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2138
2139 /* Add two lo words. */
2140 t1 = tail_a + tail_b;
2141 bv = t1 - tail_a;
2142 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2143
2144 s2 += t1;
2145
2146 /* Renormalize (s1, s2) to (t1, s2) */
2147 t1 = s1 + s2;
2148 s2 = s2 - (t1 - s1);
2149
2150 t2 += s2;
2151
2152 /* Renormalize (t1, t2) */
2153 head_t = t1 + t2;
2154 tail_t = t2 - (head_t - t1);
2155 }
2156 head_sum[0] = head_t;
2157 tail_sum[0] = tail_t;
2158 /* Imaginary part */
2159 head_a = head_sum[1];
2160 tail_a = tail_sum[1];
2161 head_b = head_prod[1];
2162 tail_b = tail_prod[1];
2163 {
2164 /* Compute double-double = double-double + double-double. */
2165 double bv;
2166 double s1, s2, t1, t2;
2167
2168 /* Add two hi words. */
2169 s1 = head_a + head_b;
2170 bv = s1 - head_a;
2171 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2172
2173 /* Add two lo words. */
2174 t1 = tail_a + tail_b;
2175 bv = t1 - tail_a;
2176 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2177
2178 s2 += t1;
2179
2180 /* Renormalize (s1, s2) to (t1, s2) */
2181 t1 = s1 + s2;
2182 s2 = s2 - (t1 - s1);
2183
2184 t2 += s2;
2185
2186 /* Renormalize (t1, t2) */
2187 head_t = t1 + t2;
2188 tail_t = t2 - (head_t - t1);
2189 }
2190 head_sum[1] = head_t;
2191 tail_sum[1] = tail_t;
2192 }
2193 x_elem = tail_x_i[jx];
2194 {
2195 head_prod[0] = (double) a_elem[0] * x_elem;
2196 tail_prod[0] = 0.0;
2197 head_prod[1] = (double) a_elem[1] * x_elem;
2198 tail_prod[1] = 0.0;
2199 }
2200 {
2201 double head_t, tail_t;
2202 double head_a, tail_a;
2203 double head_b, tail_b;
2204 /* Real part */
2205 head_a = head_sum2[0];
2206 tail_a = tail_sum2[0];
2207 head_b = head_prod[0];
2208 tail_b = tail_prod[0];
2209 {
2210 /* Compute double-double = double-double + double-double. */
2211 double bv;
2212 double s1, s2, t1, t2;
2213
2214 /* Add two hi words. */
2215 s1 = head_a + head_b;
2216 bv = s1 - head_a;
2217 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2218
2219 /* Add two lo words. */
2220 t1 = tail_a + tail_b;
2221 bv = t1 - tail_a;
2222 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2223
2224 s2 += t1;
2225
2226 /* Renormalize (s1, s2) to (t1, s2) */
2227 t1 = s1 + s2;
2228 s2 = s2 - (t1 - s1);
2229
2230 t2 += s2;
2231
2232 /* Renormalize (t1, t2) */
2233 head_t = t1 + t2;
2234 tail_t = t2 - (head_t - t1);
2235 }
2236 head_sum2[0] = head_t;
2237 tail_sum2[0] = tail_t;
2238 /* Imaginary part */
2239 head_a = head_sum2[1];
2240 tail_a = tail_sum2[1];
2241 head_b = head_prod[1];
2242 tail_b = tail_prod[1];
2243 {
2244 /* Compute double-double = double-double + double-double. */
2245 double bv;
2246 double s1, s2, t1, t2;
2247
2248 /* Add two hi words. */
2249 s1 = head_a + head_b;
2250 bv = s1 - head_a;
2251 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2252
2253 /* Add two lo words. */
2254 t1 = tail_a + tail_b;
2255 bv = t1 - tail_a;
2256 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2257
2258 s2 += t1;
2259
2260 /* Renormalize (s1, s2) to (t1, s2) */
2261 t1 = s1 + s2;
2262 s2 = s2 - (t1 - s1);
2263
2264 t2 += s2;
2265
2266 /* Renormalize (t1, t2) */
2267 head_t = t1 + t2;
2268 tail_t = t2 - (head_t - t1);
2269 }
2270 head_sum2[1] = head_t;
2271 tail_sum2[1] = tail_t;
2272 }
2273 aij += incaij;
2274 jx += incx;
2275 }
2276 {
2277 double head_t, tail_t;
2278 double head_a, tail_a;
2279 double head_b, tail_b;
2280 /* Real part */
2281 head_a = head_sum[0];
2282 tail_a = tail_sum[0];
2283 head_b = head_sum2[0];
2284 tail_b = tail_sum2[0];
2285 {
2286 /* Compute double-double = double-double + double-double. */
2287 double bv;
2288 double s1, s2, t1, t2;
2289
2290 /* Add two hi words. */
2291 s1 = head_a + head_b;
2292 bv = s1 - head_a;
2293 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2294
2295 /* Add two lo words. */
2296 t1 = tail_a + tail_b;
2297 bv = t1 - tail_a;
2298 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2299
2300 s2 += t1;
2301
2302 /* Renormalize (s1, s2) to (t1, s2) */
2303 t1 = s1 + s2;
2304 s2 = s2 - (t1 - s1);
2305
2306 t2 += s2;
2307
2308 /* Renormalize (t1, t2) */
2309 head_t = t1 + t2;
2310 tail_t = t2 - (head_t - t1);
2311 }
2312 head_sum[0] = head_t;
2313 tail_sum[0] = tail_t;
2314 /* Imaginary part */
2315 head_a = head_sum[1];
2316 tail_a = tail_sum[1];
2317 head_b = head_sum2[1];
2318 tail_b = tail_sum2[1];
2319 {
2320 /* Compute double-double = double-double + double-double. */
2321 double bv;
2322 double s1, s2, t1, t2;
2323
2324 /* Add two hi words. */
2325 s1 = head_a + head_b;
2326 bv = s1 - head_a;
2327 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2328
2329 /* Add two lo words. */
2330 t1 = tail_a + tail_b;
2331 bv = t1 - tail_a;
2332 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2333
2334 s2 += t1;
2335
2336 /* Renormalize (s1, s2) to (t1, s2) */
2337 t1 = s1 + s2;
2338 s2 = s2 - (t1 - s1);
2339
2340 t2 += s2;
2341
2342 /* Renormalize (t1, t2) */
2343 head_t = t1 + t2;
2344 tail_t = t2 - (head_t - t1);
2345 }
2346 head_sum[1] = head_t;
2347 tail_sum[1] = tail_t;
2348 }
2349 y_elem[0] = y_i[iy];
2350 y_elem[1] = y_i[iy + 1];
2351 {
2352 double head_e1, tail_e1;
2353 double d1;
2354 double d2;
2355 /* Real part */
2356 d1 = (double) y_elem[0] * beta_i[0];
2357 d2 = (double) -y_elem[1] * beta_i[1];
2358 {
2359 /* Compute double-double = double + double. */
2360 double e, t1, t2;
2361
2362 /* Knuth trick. */
2363 t1 = d1 + d2;
2364 e = t1 - d1;
2365 t2 = ((d2 - e) + (d1 - (t1 - e)));
2366
2367 /* The result is t1 + t2, after normalization. */
2368 head_e1 = t1 + t2;
2369 tail_e1 = t2 - (head_e1 - t1);
2370 }
2371 head_tmp1[0] = head_e1;
2372 tail_tmp1[0] = tail_e1;
2373 /* imaginary part */
2374 d1 = (double) y_elem[0] * beta_i[1];
2375 d2 = (double) y_elem[1] * beta_i[0];
2376 {
2377 /* Compute double-double = double + double. */
2378 double e, t1, t2;
2379
2380 /* Knuth trick. */
2381 t1 = d1 + d2;
2382 e = t1 - d1;
2383 t2 = ((d2 - e) + (d1 - (t1 - e)));
2384
2385 /* The result is t1 + t2, after normalization. */
2386 head_e1 = t1 + t2;
2387 tail_e1 = t2 - (head_e1 - t1);
2388 }
2389 head_tmp1[1] = head_e1;
2390 tail_tmp1[1] = tail_e1;
2391 }
2392 {
2393 double head_t, tail_t;
2394 double head_a, tail_a;
2395 double head_b, tail_b;
2396 /* Real part */
2397 head_a = head_sum[0];
2398 tail_a = tail_sum[0];
2399 head_b = head_tmp1[0];
2400 tail_b = tail_tmp1[0];
2401 {
2402 /* Compute double-double = double-double + double-double. */
2403 double bv;
2404 double s1, s2, t1, t2;
2405
2406 /* Add two hi words. */
2407 s1 = head_a + head_b;
2408 bv = s1 - head_a;
2409 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2410
2411 /* Add two lo words. */
2412 t1 = tail_a + tail_b;
2413 bv = t1 - tail_a;
2414 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2415
2416 s2 += t1;
2417
2418 /* Renormalize (s1, s2) to (t1, s2) */
2419 t1 = s1 + s2;
2420 s2 = s2 - (t1 - s1);
2421
2422 t2 += s2;
2423
2424 /* Renormalize (t1, t2) */
2425 head_t = t1 + t2;
2426 tail_t = t2 - (head_t - t1);
2427 }
2428 head_tmp2[0] = head_t;
2429 tail_tmp2[0] = tail_t;
2430 /* Imaginary part */
2431 head_a = head_sum[1];
2432 tail_a = tail_sum[1];
2433 head_b = head_tmp1[1];
2434 tail_b = tail_tmp1[1];
2435 {
2436 /* Compute double-double = double-double + double-double. */
2437 double bv;
2438 double s1, s2, t1, t2;
2439
2440 /* Add two hi words. */
2441 s1 = head_a + head_b;
2442 bv = s1 - head_a;
2443 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2444
2445 /* Add two lo words. */
2446 t1 = tail_a + tail_b;
2447 bv = t1 - tail_a;
2448 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2449
2450 s2 += t1;
2451
2452 /* Renormalize (s1, s2) to (t1, s2) */
2453 t1 = s1 + s2;
2454 s2 = s2 - (t1 - s1);
2455
2456 t2 += s2;
2457
2458 /* Renormalize (t1, t2) */
2459 head_t = t1 + t2;
2460 tail_t = t2 - (head_t - t1);
2461 }
2462 head_tmp2[1] = head_t;
2463 tail_tmp2[1] = tail_t;
2464 }
2465 y_i[iy] = head_tmp2[0];
2466 y_i[iy + 1] = head_tmp2[1];
2467 ai += incai;
2468 iy += incy;
2469 }
2470 } else { /* alpha != 1, the most general form:
2471 y = alpha*A*head_x + alpha*A*tail_x + beta*y */
2472 ai = 0;
2473 iy = ky;
2474 for (i = 0; i < leny; i++) {
2475 head_sum[0] = head_sum[1] = tail_sum[0] = tail_sum[1] = 0.0;;
2476 head_sum2[0] = head_sum2[1] = tail_sum2[0] = tail_sum2[1] =
2477 0.0;;
2478 aij = ai;
2479 jx = kx;
2480 for (j = 0; j < lenx; j++) {
2481 a_elem[0] = a_i[aij];
2482 a_elem[1] = a_i[aij + 1];
2483 a_elem[1] = -a_elem[1];
2484 x_elem = head_x_i[jx];
2485 {
2486 head_prod[0] = (double) a_elem[0] * x_elem;
2487 tail_prod[0] = 0.0;
2488 head_prod[1] = (double) a_elem[1] * x_elem;
2489 tail_prod[1] = 0.0;
2490 }
2491 {
2492 double head_t, tail_t;
2493 double head_a, tail_a;
2494 double head_b, tail_b;
2495 /* Real part */
2496 head_a = head_sum[0];
2497 tail_a = tail_sum[0];
2498 head_b = head_prod[0];
2499 tail_b = tail_prod[0];
2500 {
2501 /* Compute double-double = double-double + double-double. */
2502 double bv;
2503 double s1, s2, t1, t2;
2504
2505 /* Add two hi words. */
2506 s1 = head_a + head_b;
2507 bv = s1 - head_a;
2508 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2509
2510 /* Add two lo words. */
2511 t1 = tail_a + tail_b;
2512 bv = t1 - tail_a;
2513 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2514
2515 s2 += t1;
2516
2517 /* Renormalize (s1, s2) to (t1, s2) */
2518 t1 = s1 + s2;
2519 s2 = s2 - (t1 - s1);
2520
2521 t2 += s2;
2522
2523 /* Renormalize (t1, t2) */
2524 head_t = t1 + t2;
2525 tail_t = t2 - (head_t - t1);
2526 }
2527 head_sum[0] = head_t;
2528 tail_sum[0] = tail_t;
2529 /* Imaginary part */
2530 head_a = head_sum[1];
2531 tail_a = tail_sum[1];
2532 head_b = head_prod[1];
2533 tail_b = tail_prod[1];
2534 {
2535 /* Compute double-double = double-double + double-double. */
2536 double bv;
2537 double s1, s2, t1, t2;
2538
2539 /* Add two hi words. */
2540 s1 = head_a + head_b;
2541 bv = s1 - head_a;
2542 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2543
2544 /* Add two lo words. */
2545 t1 = tail_a + tail_b;
2546 bv = t1 - tail_a;
2547 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2548
2549 s2 += t1;
2550
2551 /* Renormalize (s1, s2) to (t1, s2) */
2552 t1 = s1 + s2;
2553 s2 = s2 - (t1 - s1);
2554
2555 t2 += s2;
2556
2557 /* Renormalize (t1, t2) */
2558 head_t = t1 + t2;
2559 tail_t = t2 - (head_t - t1);
2560 }
2561 head_sum[1] = head_t;
2562 tail_sum[1] = tail_t;
2563 }
2564 x_elem = tail_x_i[jx];
2565 {
2566 head_prod[0] = (double) a_elem[0] * x_elem;
2567 tail_prod[0] = 0.0;
2568 head_prod[1] = (double) a_elem[1] * x_elem;
2569 tail_prod[1] = 0.0;
2570 }
2571 {
2572 double head_t, tail_t;
2573 double head_a, tail_a;
2574 double head_b, tail_b;
2575 /* Real part */
2576 head_a = head_sum2[0];
2577 tail_a = tail_sum2[0];
2578 head_b = head_prod[0];
2579 tail_b = tail_prod[0];
2580 {
2581 /* Compute double-double = double-double + double-double. */
2582 double bv;
2583 double s1, s2, t1, t2;
2584
2585 /* Add two hi words. */
2586 s1 = head_a + head_b;
2587 bv = s1 - head_a;
2588 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2589
2590 /* Add two lo words. */
2591 t1 = tail_a + tail_b;
2592 bv = t1 - tail_a;
2593 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2594
2595 s2 += t1;
2596
2597 /* Renormalize (s1, s2) to (t1, s2) */
2598 t1 = s1 + s2;
2599 s2 = s2 - (t1 - s1);
2600
2601 t2 += s2;
2602
2603 /* Renormalize (t1, t2) */
2604 head_t = t1 + t2;
2605 tail_t = t2 - (head_t - t1);
2606 }
2607 head_sum2[0] = head_t;
2608 tail_sum2[0] = tail_t;
2609 /* Imaginary part */
2610 head_a = head_sum2[1];
2611 tail_a = tail_sum2[1];
2612 head_b = head_prod[1];
2613 tail_b = tail_prod[1];
2614 {
2615 /* Compute double-double = double-double + double-double. */
2616 double bv;
2617 double s1, s2, t1, t2;
2618
2619 /* Add two hi words. */
2620 s1 = head_a + head_b;
2621 bv = s1 - head_a;
2622 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2623
2624 /* Add two lo words. */
2625 t1 = tail_a + tail_b;
2626 bv = t1 - tail_a;
2627 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2628
2629 s2 += t1;
2630
2631 /* Renormalize (s1, s2) to (t1, s2) */
2632 t1 = s1 + s2;
2633 s2 = s2 - (t1 - s1);
2634
2635 t2 += s2;
2636
2637 /* Renormalize (t1, t2) */
2638 head_t = t1 + t2;
2639 tail_t = t2 - (head_t - t1);
2640 }
2641 head_sum2[1] = head_t;
2642 tail_sum2[1] = tail_t;
2643 }
2644 aij += incaij;
2645 jx += incx;
2646 }
2647 {
2648 double cd[2];
2649 cd[0] = (double) alpha_i[0];
2650 cd[1] = (double) alpha_i[1];
2651 {
2652 /* Compute complex-extra = complex-extra * complex-double. */
2653 double head_a0, tail_a0;
2654 double head_a1, tail_a1;
2655 double head_t1, tail_t1;
2656 double head_t2, tail_t2;
2657 head_a0 = head_sum[0];
2658 tail_a0 = tail_sum[0];
2659 head_a1 = head_sum[1];
2660 tail_a1 = tail_sum[1];
2661 /* real part */
2662 {
2663 /* Compute double-double = double-double * double. */
2664 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
2665
2666 con = head_a0 * split;
2667 a11 = con - head_a0;
2668 a11 = con - a11;
2669 a21 = head_a0 - a11;
2670 con = cd[0] * split;
2671 b1 = con - cd[0];
2672 b1 = con - b1;
2673 b2 = cd[0] - b1;
2674
2675 c11 = head_a0 * cd[0];
2676 c21 =
2677 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
2678
2679 c2 = tail_a0 * cd[0];
2680 t1 = c11 + c2;
2681 t2 = (c2 - (t1 - c11)) + c21;
2682
2683 head_t1 = t1 + t2;
2684 tail_t1 = t2 - (head_t1 - t1);
2685 }
2686 {
2687 /* Compute double-double = double-double * double. */
2688 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
2689
2690 con = head_a1 * split;
2691 a11 = con - head_a1;
2692 a11 = con - a11;
2693 a21 = head_a1 - a11;
2694 con = cd[1] * split;
2695 b1 = con - cd[1];
2696 b1 = con - b1;
2697 b2 = cd[1] - b1;
2698
2699 c11 = head_a1 * cd[1];
2700 c21 =
2701 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
2702
2703 c2 = tail_a1 * cd[1];
2704 t1 = c11 + c2;
2705 t2 = (c2 - (t1 - c11)) + c21;
2706
2707 head_t2 = t1 + t2;
2708 tail_t2 = t2 - (head_t2 - t1);
2709 }
2710 head_t2 = -head_t2;
2711 tail_t2 = -tail_t2;
2712 {
2713 /* Compute double-double = double-double + double-double. */
2714 double bv;
2715 double s1, s2, t1, t2;
2716
2717 /* Add two hi words. */
2718 s1 = head_t1 + head_t2;
2719 bv = s1 - head_t1;
2720 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
2721
2722 /* Add two lo words. */
2723 t1 = tail_t1 + tail_t2;
2724 bv = t1 - tail_t1;
2725 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
2726
2727 s2 += t1;
2728
2729 /* Renormalize (s1, s2) to (t1, s2) */
2730 t1 = s1 + s2;
2731 s2 = s2 - (t1 - s1);
2732
2733 t2 += s2;
2734
2735 /* Renormalize (t1, t2) */
2736 head_t1 = t1 + t2;
2737 tail_t1 = t2 - (head_t1 - t1);
2738 }
2739 head_tmp1[0] = head_t1;
2740 tail_tmp1[0] = tail_t1;
2741 /* imaginary part */
2742 {
2743 /* Compute double-double = double-double * double. */
2744 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
2745
2746 con = head_a1 * split;
2747 a11 = con - head_a1;
2748 a11 = con - a11;
2749 a21 = head_a1 - a11;
2750 con = cd[0] * split;
2751 b1 = con - cd[0];
2752 b1 = con - b1;
2753 b2 = cd[0] - b1;
2754
2755 c11 = head_a1 * cd[0];
2756 c21 =
2757 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
2758
2759 c2 = tail_a1 * cd[0];
2760 t1 = c11 + c2;
2761 t2 = (c2 - (t1 - c11)) + c21;
2762
2763 head_t1 = t1 + t2;
2764 tail_t1 = t2 - (head_t1 - t1);
2765 }
2766 {
2767 /* Compute double-double = double-double * double. */
2768 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
2769
2770 con = head_a0 * split;
2771 a11 = con - head_a0;
2772 a11 = con - a11;
2773 a21 = head_a0 - a11;
2774 con = cd[1] * split;
2775 b1 = con - cd[1];
2776 b1 = con - b1;
2777 b2 = cd[1] - b1;
2778
2779 c11 = head_a0 * cd[1];
2780 c21 =
2781 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
2782
2783 c2 = tail_a0 * cd[1];
2784 t1 = c11 + c2;
2785 t2 = (c2 - (t1 - c11)) + c21;
2786
2787 head_t2 = t1 + t2;
2788 tail_t2 = t2 - (head_t2 - t1);
2789 }
2790 {
2791 /* Compute double-double = double-double + double-double. */
2792 double bv;
2793 double s1, s2, t1, t2;
2794
2795 /* Add two hi words. */
2796 s1 = head_t1 + head_t2;
2797 bv = s1 - head_t1;
2798 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
2799
2800 /* Add two lo words. */
2801 t1 = tail_t1 + tail_t2;
2802 bv = t1 - tail_t1;
2803 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
2804
2805 s2 += t1;
2806
2807 /* Renormalize (s1, s2) to (t1, s2) */
2808 t1 = s1 + s2;
2809 s2 = s2 - (t1 - s1);
2810
2811 t2 += s2;
2812
2813 /* Renormalize (t1, t2) */
2814 head_t1 = t1 + t2;
2815 tail_t1 = t2 - (head_t1 - t1);
2816 }
2817 head_tmp1[1] = head_t1;
2818 tail_tmp1[1] = tail_t1;
2819 }
2820
2821 }
2822 {
2823 double cd[2];
2824 cd[0] = (double) alpha_i[0];
2825 cd[1] = (double) alpha_i[1];
2826 {
2827 /* Compute complex-extra = complex-extra * complex-double. */
2828 double head_a0, tail_a0;
2829 double head_a1, tail_a1;
2830 double head_t1, tail_t1;
2831 double head_t2, tail_t2;
2832 head_a0 = head_sum2[0];
2833 tail_a0 = tail_sum2[0];
2834 head_a1 = head_sum2[1];
2835 tail_a1 = tail_sum2[1];
2836 /* real part */
2837 {
2838 /* Compute double-double = double-double * double. */
2839 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
2840
2841 con = head_a0 * split;
2842 a11 = con - head_a0;
2843 a11 = con - a11;
2844 a21 = head_a0 - a11;
2845 con = cd[0] * split;
2846 b1 = con - cd[0];
2847 b1 = con - b1;
2848 b2 = cd[0] - b1;
2849
2850 c11 = head_a0 * cd[0];
2851 c21 =
2852 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
2853
2854 c2 = tail_a0 * cd[0];
2855 t1 = c11 + c2;
2856 t2 = (c2 - (t1 - c11)) + c21;
2857
2858 head_t1 = t1 + t2;
2859 tail_t1 = t2 - (head_t1 - t1);
2860 }
2861 {
2862 /* Compute double-double = double-double * double. */
2863 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
2864
2865 con = head_a1 * split;
2866 a11 = con - head_a1;
2867 a11 = con - a11;
2868 a21 = head_a1 - a11;
2869 con = cd[1] * split;
2870 b1 = con - cd[1];
2871 b1 = con - b1;
2872 b2 = cd[1] - b1;
2873
2874 c11 = head_a1 * cd[1];
2875 c21 =
2876 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
2877
2878 c2 = tail_a1 * cd[1];
2879 t1 = c11 + c2;
2880 t2 = (c2 - (t1 - c11)) + c21;
2881
2882 head_t2 = t1 + t2;
2883 tail_t2 = t2 - (head_t2 - t1);
2884 }
2885 head_t2 = -head_t2;
2886 tail_t2 = -tail_t2;
2887 {
2888 /* Compute double-double = double-double + double-double. */
2889 double bv;
2890 double s1, s2, t1, t2;
2891
2892 /* Add two hi words. */
2893 s1 = head_t1 + head_t2;
2894 bv = s1 - head_t1;
2895 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
2896
2897 /* Add two lo words. */
2898 t1 = tail_t1 + tail_t2;
2899 bv = t1 - tail_t1;
2900 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
2901
2902 s2 += t1;
2903
2904 /* Renormalize (s1, s2) to (t1, s2) */
2905 t1 = s1 + s2;
2906 s2 = s2 - (t1 - s1);
2907
2908 t2 += s2;
2909
2910 /* Renormalize (t1, t2) */
2911 head_t1 = t1 + t2;
2912 tail_t1 = t2 - (head_t1 - t1);
2913 }
2914 head_tmp2[0] = head_t1;
2915 tail_tmp2[0] = tail_t1;
2916 /* imaginary part */
2917 {
2918 /* Compute double-double = double-double * double. */
2919 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
2920
2921 con = head_a1 * split;
2922 a11 = con - head_a1;
2923 a11 = con - a11;
2924 a21 = head_a1 - a11;
2925 con = cd[0] * split;
2926 b1 = con - cd[0];
2927 b1 = con - b1;
2928 b2 = cd[0] - b1;
2929
2930 c11 = head_a1 * cd[0];
2931 c21 =
2932 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
2933
2934 c2 = tail_a1 * cd[0];
2935 t1 = c11 + c2;
2936 t2 = (c2 - (t1 - c11)) + c21;
2937
2938 head_t1 = t1 + t2;
2939 tail_t1 = t2 - (head_t1 - t1);
2940 }
2941 {
2942 /* Compute double-double = double-double * double. */
2943 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
2944
2945 con = head_a0 * split;
2946 a11 = con - head_a0;
2947 a11 = con - a11;
2948 a21 = head_a0 - a11;
2949 con = cd[1] * split;
2950 b1 = con - cd[1];
2951 b1 = con - b1;
2952 b2 = cd[1] - b1;
2953
2954 c11 = head_a0 * cd[1];
2955 c21 =
2956 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
2957
2958 c2 = tail_a0 * cd[1];
2959 t1 = c11 + c2;
2960 t2 = (c2 - (t1 - c11)) + c21;
2961
2962 head_t2 = t1 + t2;
2963 tail_t2 = t2 - (head_t2 - t1);
2964 }
2965 {
2966 /* Compute double-double = double-double + double-double. */
2967 double bv;
2968 double s1, s2, t1, t2;
2969
2970 /* Add two hi words. */
2971 s1 = head_t1 + head_t2;
2972 bv = s1 - head_t1;
2973 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
2974
2975 /* Add two lo words. */
2976 t1 = tail_t1 + tail_t2;
2977 bv = t1 - tail_t1;
2978 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
2979
2980 s2 += t1;
2981
2982 /* Renormalize (s1, s2) to (t1, s2) */
2983 t1 = s1 + s2;
2984 s2 = s2 - (t1 - s1);
2985
2986 t2 += s2;
2987
2988 /* Renormalize (t1, t2) */
2989 head_t1 = t1 + t2;
2990 tail_t1 = t2 - (head_t1 - t1);
2991 }
2992 head_tmp2[1] = head_t1;
2993 tail_tmp2[1] = tail_t1;
2994 }
2995
2996 }
2997 {
2998 double head_t, tail_t;
2999 double head_a, tail_a;
3000 double head_b, tail_b;
3001 /* Real part */
3002 head_a = head_tmp1[0];
3003 tail_a = tail_tmp1[0];
3004 head_b = head_tmp2[0];
3005 tail_b = tail_tmp2[0];
3006 {
3007 /* Compute double-double = double-double + double-double. */
3008 double bv;
3009 double s1, s2, t1, t2;
3010
3011 /* Add two hi words. */
3012 s1 = head_a + head_b;
3013 bv = s1 - head_a;
3014 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
3015
3016 /* Add two lo words. */
3017 t1 = tail_a + tail_b;
3018 bv = t1 - tail_a;
3019 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
3020
3021 s2 += t1;
3022
3023 /* Renormalize (s1, s2) to (t1, s2) */
3024 t1 = s1 + s2;
3025 s2 = s2 - (t1 - s1);
3026
3027 t2 += s2;
3028
3029 /* Renormalize (t1, t2) */
3030 head_t = t1 + t2;
3031 tail_t = t2 - (head_t - t1);
3032 }
3033 head_tmp1[0] = head_t;
3034 tail_tmp1[0] = tail_t;
3035 /* Imaginary part */
3036 head_a = head_tmp1[1];
3037 tail_a = tail_tmp1[1];
3038 head_b = head_tmp2[1];
3039 tail_b = tail_tmp2[1];
3040 {
3041 /* Compute double-double = double-double + double-double. */
3042 double bv;
3043 double s1, s2, t1, t2;
3044
3045 /* Add two hi words. */
3046 s1 = head_a + head_b;
3047 bv = s1 - head_a;
3048 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
3049
3050 /* Add two lo words. */
3051 t1 = tail_a + tail_b;
3052 bv = t1 - tail_a;
3053 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
3054
3055 s2 += t1;
3056
3057 /* Renormalize (s1, s2) to (t1, s2) */
3058 t1 = s1 + s2;
3059 s2 = s2 - (t1 - s1);
3060
3061 t2 += s2;
3062
3063 /* Renormalize (t1, t2) */
3064 head_t = t1 + t2;
3065 tail_t = t2 - (head_t - t1);
3066 }
3067 head_tmp1[1] = head_t;
3068 tail_tmp1[1] = tail_t;
3069 }
3070 y_elem[0] = y_i[iy];
3071 y_elem[1] = y_i[iy + 1];
3072 {
3073 double head_e1, tail_e1;
3074 double d1;
3075 double d2;
3076 /* Real part */
3077 d1 = (double) y_elem[0] * beta_i[0];
3078 d2 = (double) -y_elem[1] * beta_i[1];
3079 {
3080 /* Compute double-double = double + double. */
3081 double e, t1, t2;
3082
3083 /* Knuth trick. */
3084 t1 = d1 + d2;
3085 e = t1 - d1;
3086 t2 = ((d2 - e) + (d1 - (t1 - e)));
3087
3088 /* The result is t1 + t2, after normalization. */
3089 head_e1 = t1 + t2;
3090 tail_e1 = t2 - (head_e1 - t1);
3091 }
3092 head_tmp2[0] = head_e1;
3093 tail_tmp2[0] = tail_e1;
3094 /* imaginary part */
3095 d1 = (double) y_elem[0] * beta_i[1];
3096 d2 = (double) y_elem[1] * beta_i[0];
3097 {
3098 /* Compute double-double = double + double. */
3099 double e, t1, t2;
3100
3101 /* Knuth trick. */
3102 t1 = d1 + d2;
3103 e = t1 - d1;
3104 t2 = ((d2 - e) + (d1 - (t1 - e)));
3105
3106 /* The result is t1 + t2, after normalization. */
3107 head_e1 = t1 + t2;
3108 tail_e1 = t2 - (head_e1 - t1);
3109 }
3110 head_tmp2[1] = head_e1;
3111 tail_tmp2[1] = tail_e1;
3112 }
3113 {
3114 double head_t, tail_t;
3115 double head_a, tail_a;
3116 double head_b, tail_b;
3117 /* Real part */
3118 head_a = head_tmp1[0];
3119 tail_a = tail_tmp1[0];
3120 head_b = head_tmp2[0];
3121 tail_b = tail_tmp2[0];
3122 {
3123 /* Compute double-double = double-double + double-double. */
3124 double bv;
3125 double s1, s2, t1, t2;
3126
3127 /* Add two hi words. */
3128 s1 = head_a + head_b;
3129 bv = s1 - head_a;
3130 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
3131
3132 /* Add two lo words. */
3133 t1 = tail_a + tail_b;
3134 bv = t1 - tail_a;
3135 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
3136
3137 s2 += t1;
3138
3139 /* Renormalize (s1, s2) to (t1, s2) */
3140 t1 = s1 + s2;
3141 s2 = s2 - (t1 - s1);
3142
3143 t2 += s2;
3144
3145 /* Renormalize (t1, t2) */
3146 head_t = t1 + t2;
3147 tail_t = t2 - (head_t - t1);
3148 }
3149 head_tmp1[0] = head_t;
3150 tail_tmp1[0] = tail_t;
3151 /* Imaginary part */
3152 head_a = head_tmp1[1];
3153 tail_a = tail_tmp1[1];
3154 head_b = head_tmp2[1];
3155 tail_b = tail_tmp2[1];
3156 {
3157 /* Compute double-double = double-double + double-double. */
3158 double bv;
3159 double s1, s2, t1, t2;
3160
3161 /* Add two hi words. */
3162 s1 = head_a + head_b;
3163 bv = s1 - head_a;
3164 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
3165
3166 /* Add two lo words. */
3167 t1 = tail_a + tail_b;
3168 bv = t1 - tail_a;
3169 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
3170
3171 s2 += t1;
3172
3173 /* Renormalize (s1, s2) to (t1, s2) */
3174 t1 = s1 + s2;
3175 s2 = s2 - (t1 - s1);
3176
3177 t2 += s2;
3178
3179 /* Renormalize (t1, t2) */
3180 head_t = t1 + t2;
3181 tail_t = t2 - (head_t - t1);
3182 }
3183 head_tmp1[1] = head_t;
3184 tail_tmp1[1] = tail_t;
3185 }
3186 y_i[iy] = head_tmp1[0];
3187 y_i[iy + 1] = head_tmp1[1];
3188 ai += incai;
3189 iy += incy;
3190 }
3191 }
3192 }
3193
3194 } else {
3195
3196 /* if beta = 0, we can save m multiplies:
3197 y = alpha*A*head_x + alpha*A*tail_x */
3198 if (beta_i[0] == 0.0 && beta_i[1] == 0.0) {
3199 if ((alpha_i[0] == 1.0 && alpha_i[1] == 0.0)) {
3200 /* save m more multiplies if alpha = 1 */
3201 ai = 0;
3202 iy = ky;
3203 for (i = 0; i < leny; i++) {
3204 head_sum[0] = head_sum[1] = tail_sum[0] = tail_sum[1] = 0.0;
3205 head_sum2[0] = head_sum2[1] = tail_sum2[0] = tail_sum2[1] =
3206 0.0;
3207 aij = ai;
3208 jx = kx;
3209 for (j = 0; j < lenx; j++) {
3210 a_elem[0] = a_i[aij];
3211 a_elem[1] = a_i[aij + 1];
3212
3213 x_elem = head_x_i[jx];
3214 {
3215 head_prod[0] = (double) a_elem[0] * x_elem;
3216 tail_prod[0] = 0.0;
3217 head_prod[1] = (double) a_elem[1] * x_elem;
3218 tail_prod[1] = 0.0;
3219 }
3220 {
3221 double head_t, tail_t;
3222 double head_a, tail_a;
3223 double head_b, tail_b;
3224 /* Real part */
3225 head_a = head_sum[0];
3226 tail_a = tail_sum[0];
3227 head_b = head_prod[0];
3228 tail_b = tail_prod[0];
3229 {
3230 /* Compute double-double = double-double + double-double. */
3231 double bv;
3232 double s1, s2, t1, t2;
3233
3234 /* Add two hi words. */
3235 s1 = head_a + head_b;
3236 bv = s1 - head_a;
3237 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
3238
3239 /* Add two lo words. */
3240 t1 = tail_a + tail_b;
3241 bv = t1 - tail_a;
3242 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
3243
3244 s2 += t1;
3245
3246 /* Renormalize (s1, s2) to (t1, s2) */
3247 t1 = s1 + s2;
3248 s2 = s2 - (t1 - s1);
3249
3250 t2 += s2;
3251
3252 /* Renormalize (t1, t2) */
3253 head_t = t1 + t2;
3254 tail_t = t2 - (head_t - t1);
3255 }
3256 head_sum[0] = head_t;
3257 tail_sum[0] = tail_t;
3258 /* Imaginary part */
3259 head_a = head_sum[1];
3260 tail_a = tail_sum[1];
3261 head_b = head_prod[1];
3262 tail_b = tail_prod[1];
3263 {
3264 /* Compute double-double = double-double + double-double. */
3265 double bv;
3266 double s1, s2, t1, t2;
3267
3268 /* Add two hi words. */
3269 s1 = head_a + head_b;
3270 bv = s1 - head_a;
3271 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
3272
3273 /* Add two lo words. */
3274 t1 = tail_a + tail_b;
3275 bv = t1 - tail_a;
3276 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
3277
3278 s2 += t1;
3279
3280 /* Renormalize (s1, s2) to (t1, s2) */
3281 t1 = s1 + s2;
3282 s2 = s2 - (t1 - s1);
3283
3284 t2 += s2;
3285
3286 /* Renormalize (t1, t2) */
3287 head_t = t1 + t2;
3288 tail_t = t2 - (head_t - t1);
3289 }
3290 head_sum[1] = head_t;
3291 tail_sum[1] = tail_t;
3292 }
3293 x_elem = tail_x_i[jx];
3294 {
3295 head_prod[0] = (double) a_elem[0] * x_elem;
3296 tail_prod[0] = 0.0;
3297 head_prod[1] = (double) a_elem[1] * x_elem;
3298 tail_prod[1] = 0.0;
3299 }
3300 {
3301 double head_t, tail_t;
3302 double head_a, tail_a;
3303 double head_b, tail_b;
3304 /* Real part */
3305 head_a = head_sum2[0];
3306 tail_a = tail_sum2[0];
3307 head_b = head_prod[0];
3308 tail_b = tail_prod[0];
3309 {
3310 /* Compute double-double = double-double + double-double. */
3311 double bv;
3312 double s1, s2, t1, t2;
3313
3314 /* Add two hi words. */
3315 s1 = head_a + head_b;
3316 bv = s1 - head_a;
3317 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
3318
3319 /* Add two lo words. */
3320 t1 = tail_a + tail_b;
3321 bv = t1 - tail_a;
3322 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
3323
3324 s2 += t1;
3325
3326 /* Renormalize (s1, s2) to (t1, s2) */
3327 t1 = s1 + s2;
3328 s2 = s2 - (t1 - s1);
3329
3330 t2 += s2;
3331
3332 /* Renormalize (t1, t2) */
3333 head_t = t1 + t2;
3334 tail_t = t2 - (head_t - t1);
3335 }
3336 head_sum2[0] = head_t;
3337 tail_sum2[0] = tail_t;
3338 /* Imaginary part */
3339 head_a = head_sum2[1];
3340 tail_a = tail_sum2[1];
3341 head_b = head_prod[1];
3342 tail_b = tail_prod[1];
3343 {
3344 /* Compute double-double = double-double + double-double. */
3345 double bv;
3346 double s1, s2, t1, t2;
3347
3348 /* Add two hi words. */
3349 s1 = head_a + head_b;
3350 bv = s1 - head_a;
3351 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
3352
3353 /* Add two lo words. */
3354 t1 = tail_a + tail_b;
3355 bv = t1 - tail_a;
3356 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
3357
3358 s2 += t1;
3359
3360 /* Renormalize (s1, s2) to (t1, s2) */
3361 t1 = s1 + s2;
3362 s2 = s2 - (t1 - s1);
3363
3364 t2 += s2;
3365
3366 /* Renormalize (t1, t2) */
3367 head_t = t1 + t2;
3368 tail_t = t2 - (head_t - t1);
3369 }
3370 head_sum2[1] = head_t;
3371 tail_sum2[1] = tail_t;
3372 }
3373 aij += incaij;
3374 jx += incx;
3375 }
3376 {
3377 double head_t, tail_t;
3378 double head_a, tail_a;
3379 double head_b, tail_b;
3380 /* Real part */
3381 head_a = head_sum[0];
3382 tail_a = tail_sum[0];
3383 head_b = head_sum2[0];
3384 tail_b = tail_sum2[0];
3385 {
3386 /* Compute double-double = double-double + double-double. */
3387 double bv;
3388 double s1, s2, t1, t2;
3389
3390 /* Add two hi words. */
3391 s1 = head_a + head_b;
3392 bv = s1 - head_a;
3393 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
3394
3395 /* Add two lo words. */
3396 t1 = tail_a + tail_b;
3397 bv = t1 - tail_a;
3398 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
3399
3400 s2 += t1;
3401
3402 /* Renormalize (s1, s2) to (t1, s2) */
3403 t1 = s1 + s2;
3404 s2 = s2 - (t1 - s1);
3405
3406 t2 += s2;
3407
3408 /* Renormalize (t1, t2) */
3409 head_t = t1 + t2;
3410 tail_t = t2 - (head_t - t1);
3411 }
3412 head_sum[0] = head_t;
3413 tail_sum[0] = tail_t;
3414 /* Imaginary part */
3415 head_a = head_sum[1];
3416 tail_a = tail_sum[1];
3417 head_b = head_sum2[1];
3418 tail_b = tail_sum2[1];
3419 {
3420 /* Compute double-double = double-double + double-double. */
3421 double bv;
3422 double s1, s2, t1, t2;
3423
3424 /* Add two hi words. */
3425 s1 = head_a + head_b;
3426 bv = s1 - head_a;
3427 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
3428
3429 /* Add two lo words. */
3430 t1 = tail_a + tail_b;
3431 bv = t1 - tail_a;
3432 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
3433
3434 s2 += t1;
3435
3436 /* Renormalize (s1, s2) to (t1, s2) */
3437 t1 = s1 + s2;
3438 s2 = s2 - (t1 - s1);
3439
3440 t2 += s2;
3441
3442 /* Renormalize (t1, t2) */
3443 head_t = t1 + t2;
3444 tail_t = t2 - (head_t - t1);
3445 }
3446 head_sum[1] = head_t;
3447 tail_sum[1] = tail_t;
3448 }
3449 y_i[iy] = head_sum[0];
3450 y_i[iy + 1] = head_sum[1];
3451 ai += incai;
3452 iy += incy;
3453 } /* end for */
3454 } else { /* alpha != 1 */
3455 ai = 0;
3456 iy = ky;
3457 for (i = 0; i < leny; i++) {
3458 head_sum[0] = head_sum[1] = tail_sum[0] = tail_sum[1] = 0.0;
3459 head_sum2[0] = head_sum2[1] = tail_sum2[0] = tail_sum2[1] =
3460 0.0;
3461 aij = ai;
3462 jx = kx;
3463 for (j = 0; j < lenx; j++) {
3464 a_elem[0] = a_i[aij];
3465 a_elem[1] = a_i[aij + 1];
3466
3467 x_elem = head_x_i[jx];
3468 {
3469 head_prod[0] = (double) a_elem[0] * x_elem;
3470 tail_prod[0] = 0.0;
3471 head_prod[1] = (double) a_elem[1] * x_elem;
3472 tail_prod[1] = 0.0;
3473 }
3474 {
3475 double head_t, tail_t;
3476 double head_a, tail_a;
3477 double head_b, tail_b;
3478 /* Real part */
3479 head_a = head_sum[0];
3480 tail_a = tail_sum[0];
3481 head_b = head_prod[0];
3482 tail_b = tail_prod[0];
3483 {
3484 /* Compute double-double = double-double + double-double. */
3485 double bv;
3486 double s1, s2, t1, t2;
3487
3488 /* Add two hi words. */
3489 s1 = head_a + head_b;
3490 bv = s1 - head_a;
3491 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
3492
3493 /* Add two lo words. */
3494 t1 = tail_a + tail_b;
3495 bv = t1 - tail_a;
3496 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
3497
3498 s2 += t1;
3499
3500 /* Renormalize (s1, s2) to (t1, s2) */
3501 t1 = s1 + s2;
3502 s2 = s2 - (t1 - s1);
3503
3504 t2 += s2;
3505
3506 /* Renormalize (t1, t2) */
3507 head_t = t1 + t2;
3508 tail_t = t2 - (head_t - t1);
3509 }
3510 head_sum[0] = head_t;
3511 tail_sum[0] = tail_t;
3512 /* Imaginary part */
3513 head_a = head_sum[1];
3514 tail_a = tail_sum[1];
3515 head_b = head_prod[1];
3516 tail_b = tail_prod[1];
3517 {
3518 /* Compute double-double = double-double + double-double. */
3519 double bv;
3520 double s1, s2, t1, t2;
3521
3522 /* Add two hi words. */
3523 s1 = head_a + head_b;
3524 bv = s1 - head_a;
3525 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
3526
3527 /* Add two lo words. */
3528 t1 = tail_a + tail_b;
3529 bv = t1 - tail_a;
3530 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
3531
3532 s2 += t1;
3533
3534 /* Renormalize (s1, s2) to (t1, s2) */
3535 t1 = s1 + s2;
3536 s2 = s2 - (t1 - s1);
3537
3538 t2 += s2;
3539
3540 /* Renormalize (t1, t2) */
3541 head_t = t1 + t2;
3542 tail_t = t2 - (head_t - t1);
3543 }
3544 head_sum[1] = head_t;
3545 tail_sum[1] = tail_t;
3546 }
3547 x_elem = tail_x_i[jx];
3548 {
3549 head_prod[0] = (double) a_elem[0] * x_elem;
3550 tail_prod[0] = 0.0;
3551 head_prod[1] = (double) a_elem[1] * x_elem;
3552 tail_prod[1] = 0.0;
3553 }
3554 {
3555 double head_t, tail_t;
3556 double head_a, tail_a;
3557 double head_b, tail_b;
3558 /* Real part */
3559 head_a = head_sum2[0];
3560 tail_a = tail_sum2[0];
3561 head_b = head_prod[0];
3562 tail_b = tail_prod[0];
3563 {
3564 /* Compute double-double = double-double + double-double. */
3565 double bv;
3566 double s1, s2, t1, t2;
3567
3568 /* Add two hi words. */
3569 s1 = head_a + head_b;
3570 bv = s1 - head_a;
3571 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
3572
3573 /* Add two lo words. */
3574 t1 = tail_a + tail_b;
3575 bv = t1 - tail_a;
3576 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
3577
3578 s2 += t1;
3579
3580 /* Renormalize (s1, s2) to (t1, s2) */
3581 t1 = s1 + s2;
3582 s2 = s2 - (t1 - s1);
3583
3584 t2 += s2;
3585
3586 /* Renormalize (t1, t2) */
3587 head_t = t1 + t2;
3588 tail_t = t2 - (head_t - t1);
3589 }
3590 head_sum2[0] = head_t;
3591 tail_sum2[0] = tail_t;
3592 /* Imaginary part */
3593 head_a = head_sum2[1];
3594 tail_a = tail_sum2[1];
3595 head_b = head_prod[1];
3596 tail_b = tail_prod[1];
3597 {
3598 /* Compute double-double = double-double + double-double. */
3599 double bv;
3600 double s1, s2, t1, t2;
3601
3602 /* Add two hi words. */
3603 s1 = head_a + head_b;
3604 bv = s1 - head_a;
3605 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
3606
3607 /* Add two lo words. */
3608 t1 = tail_a + tail_b;
3609 bv = t1 - tail_a;
3610 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
3611
3612 s2 += t1;
3613
3614 /* Renormalize (s1, s2) to (t1, s2) */
3615 t1 = s1 + s2;
3616 s2 = s2 - (t1 - s1);
3617
3618 t2 += s2;
3619
3620 /* Renormalize (t1, t2) */
3621 head_t = t1 + t2;
3622 tail_t = t2 - (head_t - t1);
3623 }
3624 head_sum2[1] = head_t;
3625 tail_sum2[1] = tail_t;
3626 }
3627 aij += incaij;
3628 jx += incx;
3629 }
3630 {
3631 double cd[2];
3632 cd[0] = (double) alpha_i[0];
3633 cd[1] = (double) alpha_i[1];
3634 {
3635 /* Compute complex-extra = complex-extra * complex-double. */
3636 double head_a0, tail_a0;
3637 double head_a1, tail_a1;
3638 double head_t1, tail_t1;
3639 double head_t2, tail_t2;
3640 head_a0 = head_sum[0];
3641 tail_a0 = tail_sum[0];
3642 head_a1 = head_sum[1];
3643 tail_a1 = tail_sum[1];
3644 /* real part */
3645 {
3646 /* Compute double-double = double-double * double. */
3647 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
3648
3649 con = head_a0 * split;
3650 a11 = con - head_a0;
3651 a11 = con - a11;
3652 a21 = head_a0 - a11;
3653 con = cd[0] * split;
3654 b1 = con - cd[0];
3655 b1 = con - b1;
3656 b2 = cd[0] - b1;
3657
3658 c11 = head_a0 * cd[0];
3659 c21 =
3660 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
3661
3662 c2 = tail_a0 * cd[0];
3663 t1 = c11 + c2;
3664 t2 = (c2 - (t1 - c11)) + c21;
3665
3666 head_t1 = t1 + t2;
3667 tail_t1 = t2 - (head_t1 - t1);
3668 }
3669 {
3670 /* Compute double-double = double-double * double. */
3671 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
3672
3673 con = head_a1 * split;
3674 a11 = con - head_a1;
3675 a11 = con - a11;
3676 a21 = head_a1 - a11;
3677 con = cd[1] * split;
3678 b1 = con - cd[1];
3679 b1 = con - b1;
3680 b2 = cd[1] - b1;
3681
3682 c11 = head_a1 * cd[1];
3683 c21 =
3684 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
3685
3686 c2 = tail_a1 * cd[1];
3687 t1 = c11 + c2;
3688 t2 = (c2 - (t1 - c11)) + c21;
3689
3690 head_t2 = t1 + t2;
3691 tail_t2 = t2 - (head_t2 - t1);
3692 }
3693 head_t2 = -head_t2;
3694 tail_t2 = -tail_t2;
3695 {
3696 /* Compute double-double = double-double + double-double. */
3697 double bv;
3698 double s1, s2, t1, t2;
3699
3700 /* Add two hi words. */
3701 s1 = head_t1 + head_t2;
3702 bv = s1 - head_t1;
3703 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
3704
3705 /* Add two lo words. */
3706 t1 = tail_t1 + tail_t2;
3707 bv = t1 - tail_t1;
3708 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
3709
3710 s2 += t1;
3711
3712 /* Renormalize (s1, s2) to (t1, s2) */
3713 t1 = s1 + s2;
3714 s2 = s2 - (t1 - s1);
3715
3716 t2 += s2;
3717
3718 /* Renormalize (t1, t2) */
3719 head_t1 = t1 + t2;
3720 tail_t1 = t2 - (head_t1 - t1);
3721 }
3722 head_tmp1[0] = head_t1;
3723 tail_tmp1[0] = tail_t1;
3724 /* imaginary part */
3725 {
3726 /* Compute double-double = double-double * double. */
3727 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
3728
3729 con = head_a1 * split;
3730 a11 = con - head_a1;
3731 a11 = con - a11;
3732 a21 = head_a1 - a11;
3733 con = cd[0] * split;
3734 b1 = con - cd[0];
3735 b1 = con - b1;
3736 b2 = cd[0] - b1;
3737
3738 c11 = head_a1 * cd[0];
3739 c21 =
3740 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
3741
3742 c2 = tail_a1 * cd[0];
3743 t1 = c11 + c2;
3744 t2 = (c2 - (t1 - c11)) + c21;
3745
3746 head_t1 = t1 + t2;
3747 tail_t1 = t2 - (head_t1 - t1);
3748 }
3749 {
3750 /* Compute double-double = double-double * double. */
3751 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
3752
3753 con = head_a0 * split;
3754 a11 = con - head_a0;
3755 a11 = con - a11;
3756 a21 = head_a0 - a11;
3757 con = cd[1] * split;
3758 b1 = con - cd[1];
3759 b1 = con - b1;
3760 b2 = cd[1] - b1;
3761
3762 c11 = head_a0 * cd[1];
3763 c21 =
3764 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
3765
3766 c2 = tail_a0 * cd[1];
3767 t1 = c11 + c2;
3768 t2 = (c2 - (t1 - c11)) + c21;
3769
3770 head_t2 = t1 + t2;
3771 tail_t2 = t2 - (head_t2 - t1);
3772 }
3773 {
3774 /* Compute double-double = double-double + double-double. */
3775 double bv;
3776 double s1, s2, t1, t2;
3777
3778 /* Add two hi words. */
3779 s1 = head_t1 + head_t2;
3780 bv = s1 - head_t1;
3781 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
3782
3783 /* Add two lo words. */
3784 t1 = tail_t1 + tail_t2;
3785 bv = t1 - tail_t1;
3786 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
3787
3788 s2 += t1;
3789
3790 /* Renormalize (s1, s2) to (t1, s2) */
3791 t1 = s1 + s2;
3792 s2 = s2 - (t1 - s1);
3793
3794 t2 += s2;
3795
3796 /* Renormalize (t1, t2) */
3797 head_t1 = t1 + t2;
3798 tail_t1 = t2 - (head_t1 - t1);
3799 }
3800 head_tmp1[1] = head_t1;
3801 tail_tmp1[1] = tail_t1;
3802 }
3803
3804 }
3805 {
3806 double cd[2];
3807 cd[0] = (double) alpha_i[0];
3808 cd[1] = (double) alpha_i[1];
3809 {
3810 /* Compute complex-extra = complex-extra * complex-double. */
3811 double head_a0, tail_a0;
3812 double head_a1, tail_a1;
3813 double head_t1, tail_t1;
3814 double head_t2, tail_t2;
3815 head_a0 = head_sum2[0];
3816 tail_a0 = tail_sum2[0];
3817 head_a1 = head_sum2[1];
3818 tail_a1 = tail_sum2[1];
3819 /* real part */
3820 {
3821 /* Compute double-double = double-double * double. */
3822 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
3823
3824 con = head_a0 * split;
3825 a11 = con - head_a0;
3826 a11 = con - a11;
3827 a21 = head_a0 - a11;
3828 con = cd[0] * split;
3829 b1 = con - cd[0];
3830 b1 = con - b1;
3831 b2 = cd[0] - b1;
3832
3833 c11 = head_a0 * cd[0];
3834 c21 =
3835 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
3836
3837 c2 = tail_a0 * cd[0];
3838 t1 = c11 + c2;
3839 t2 = (c2 - (t1 - c11)) + c21;
3840
3841 head_t1 = t1 + t2;
3842 tail_t1 = t2 - (head_t1 - t1);
3843 }
3844 {
3845 /* Compute double-double = double-double * double. */
3846 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
3847
3848 con = head_a1 * split;
3849 a11 = con - head_a1;
3850 a11 = con - a11;
3851 a21 = head_a1 - a11;
3852 con = cd[1] * split;
3853 b1 = con - cd[1];
3854 b1 = con - b1;
3855 b2 = cd[1] - b1;
3856
3857 c11 = head_a1 * cd[1];
3858 c21 =
3859 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
3860
3861 c2 = tail_a1 * cd[1];
3862 t1 = c11 + c2;
3863 t2 = (c2 - (t1 - c11)) + c21;
3864
3865 head_t2 = t1 + t2;
3866 tail_t2 = t2 - (head_t2 - t1);
3867 }
3868 head_t2 = -head_t2;
3869 tail_t2 = -tail_t2;
3870 {
3871 /* Compute double-double = double-double + double-double. */
3872 double bv;
3873 double s1, s2, t1, t2;
3874
3875 /* Add two hi words. */
3876 s1 = head_t1 + head_t2;
3877 bv = s1 - head_t1;
3878 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
3879
3880 /* Add two lo words. */
3881 t1 = tail_t1 + tail_t2;
3882 bv = t1 - tail_t1;
3883 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
3884
3885 s2 += t1;
3886
3887 /* Renormalize (s1, s2) to (t1, s2) */
3888 t1 = s1 + s2;
3889 s2 = s2 - (t1 - s1);
3890
3891 t2 += s2;
3892
3893 /* Renormalize (t1, t2) */
3894 head_t1 = t1 + t2;
3895 tail_t1 = t2 - (head_t1 - t1);
3896 }
3897 head_tmp2[0] = head_t1;
3898 tail_tmp2[0] = tail_t1;
3899 /* imaginary part */
3900 {
3901 /* Compute double-double = double-double * double. */
3902 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
3903
3904 con = head_a1 * split;
3905 a11 = con - head_a1;
3906 a11 = con - a11;
3907 a21 = head_a1 - a11;
3908 con = cd[0] * split;
3909 b1 = con - cd[0];
3910 b1 = con - b1;
3911 b2 = cd[0] - b1;
3912
3913 c11 = head_a1 * cd[0];
3914 c21 =
3915 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
3916
3917 c2 = tail_a1 * cd[0];
3918 t1 = c11 + c2;
3919 t2 = (c2 - (t1 - c11)) + c21;
3920
3921 head_t1 = t1 + t2;
3922 tail_t1 = t2 - (head_t1 - t1);
3923 }
3924 {
3925 /* Compute double-double = double-double * double. */
3926 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
3927
3928 con = head_a0 * split;
3929 a11 = con - head_a0;
3930 a11 = con - a11;
3931 a21 = head_a0 - a11;
3932 con = cd[1] * split;
3933 b1 = con - cd[1];
3934 b1 = con - b1;
3935 b2 = cd[1] - b1;
3936
3937 c11 = head_a0 * cd[1];
3938 c21 =
3939 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
3940
3941 c2 = tail_a0 * cd[1];
3942 t1 = c11 + c2;
3943 t2 = (c2 - (t1 - c11)) + c21;
3944
3945 head_t2 = t1 + t2;
3946 tail_t2 = t2 - (head_t2 - t1);
3947 }
3948 {
3949 /* Compute double-double = double-double + double-double. */
3950 double bv;
3951 double s1, s2, t1, t2;
3952
3953 /* Add two hi words. */
3954 s1 = head_t1 + head_t2;
3955 bv = s1 - head_t1;
3956 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
3957
3958 /* Add two lo words. */
3959 t1 = tail_t1 + tail_t2;
3960 bv = t1 - tail_t1;
3961 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
3962
3963 s2 += t1;
3964
3965 /* Renormalize (s1, s2) to (t1, s2) */
3966 t1 = s1 + s2;
3967 s2 = s2 - (t1 - s1);
3968
3969 t2 += s2;
3970
3971 /* Renormalize (t1, t2) */
3972 head_t1 = t1 + t2;
3973 tail_t1 = t2 - (head_t1 - t1);
3974 }
3975 head_tmp2[1] = head_t1;
3976 tail_tmp2[1] = tail_t1;
3977 }
3978
3979 }
3980 {
3981 double head_t, tail_t;
3982 double head_a, tail_a;
3983 double head_b, tail_b;
3984 /* Real part */
3985 head_a = head_tmp1[0];
3986 tail_a = tail_tmp1[0];
3987 head_b = head_tmp2[0];
3988 tail_b = tail_tmp2[0];
3989 {
3990 /* Compute double-double = double-double + double-double. */
3991 double bv;
3992 double s1, s2, t1, t2;
3993
3994 /* Add two hi words. */
3995 s1 = head_a + head_b;
3996 bv = s1 - head_a;
3997 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
3998
3999 /* Add two lo words. */
4000 t1 = tail_a + tail_b;
4001 bv = t1 - tail_a;
4002 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
4003
4004 s2 += t1;
4005
4006 /* Renormalize (s1, s2) to (t1, s2) */
4007 t1 = s1 + s2;
4008 s2 = s2 - (t1 - s1);
4009
4010 t2 += s2;
4011
4012 /* Renormalize (t1, t2) */
4013 head_t = t1 + t2;
4014 tail_t = t2 - (head_t - t1);
4015 }
4016 head_tmp1[0] = head_t;
4017 tail_tmp1[0] = tail_t;
4018 /* Imaginary part */
4019 head_a = head_tmp1[1];
4020 tail_a = tail_tmp1[1];
4021 head_b = head_tmp2[1];
4022 tail_b = tail_tmp2[1];
4023 {
4024 /* Compute double-double = double-double + double-double. */
4025 double bv;
4026 double s1, s2, t1, t2;
4027
4028 /* Add two hi words. */
4029 s1 = head_a + head_b;
4030 bv = s1 - head_a;
4031 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
4032
4033 /* Add two lo words. */
4034 t1 = tail_a + tail_b;
4035 bv = t1 - tail_a;
4036 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
4037
4038 s2 += t1;
4039
4040 /* Renormalize (s1, s2) to (t1, s2) */
4041 t1 = s1 + s2;
4042 s2 = s2 - (t1 - s1);
4043
4044 t2 += s2;
4045
4046 /* Renormalize (t1, t2) */
4047 head_t = t1 + t2;
4048 tail_t = t2 - (head_t - t1);
4049 }
4050 head_tmp1[1] = head_t;
4051 tail_tmp1[1] = tail_t;
4052 }
4053 y_i[iy] = head_tmp1[0];
4054 y_i[iy + 1] = head_tmp1[1];
4055 ai += incai;
4056 iy += incy;
4057 }
4058 }
4059 } else { /* beta != 0 */
4060 if ((alpha_i[0] == 1.0 && alpha_i[1] == 0.0)) {
4061 /* save m multiplies if alpha = 1 */
4062 ai = 0;
4063 iy = ky;
4064 for (i = 0; i < leny; i++) {
4065 head_sum[0] = head_sum[1] = tail_sum[0] = tail_sum[1] = 0.0;;
4066 head_sum2[0] = head_sum2[1] = tail_sum2[0] = tail_sum2[1] =
4067 0.0;;
4068 aij = ai;
4069 jx = kx;
4070 for (j = 0; j < lenx; j++) {
4071 a_elem[0] = a_i[aij];
4072 a_elem[1] = a_i[aij + 1];
4073
4074 x_elem = head_x_i[jx];
4075 {
4076 head_prod[0] = (double) a_elem[0] * x_elem;
4077 tail_prod[0] = 0.0;
4078 head_prod[1] = (double) a_elem[1] * x_elem;
4079 tail_prod[1] = 0.0;
4080 }
4081 {
4082 double head_t, tail_t;
4083 double head_a, tail_a;
4084 double head_b, tail_b;
4085 /* Real part */
4086 head_a = head_sum[0];
4087 tail_a = tail_sum[0];
4088 head_b = head_prod[0];
4089 tail_b = tail_prod[0];
4090 {
4091 /* Compute double-double = double-double + double-double. */
4092 double bv;
4093 double s1, s2, t1, t2;
4094
4095 /* Add two hi words. */
4096 s1 = head_a + head_b;
4097 bv = s1 - head_a;
4098 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
4099
4100 /* Add two lo words. */
4101 t1 = tail_a + tail_b;
4102 bv = t1 - tail_a;
4103 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
4104
4105 s2 += t1;
4106
4107 /* Renormalize (s1, s2) to (t1, s2) */
4108 t1 = s1 + s2;
4109 s2 = s2 - (t1 - s1);
4110
4111 t2 += s2;
4112
4113 /* Renormalize (t1, t2) */
4114 head_t = t1 + t2;
4115 tail_t = t2 - (head_t - t1);
4116 }
4117 head_sum[0] = head_t;
4118 tail_sum[0] = tail_t;
4119 /* Imaginary part */
4120 head_a = head_sum[1];
4121 tail_a = tail_sum[1];
4122 head_b = head_prod[1];
4123 tail_b = tail_prod[1];
4124 {
4125 /* Compute double-double = double-double + double-double. */
4126 double bv;
4127 double s1, s2, t1, t2;
4128
4129 /* Add two hi words. */
4130 s1 = head_a + head_b;
4131 bv = s1 - head_a;
4132 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
4133
4134 /* Add two lo words. */
4135 t1 = tail_a + tail_b;
4136 bv = t1 - tail_a;
4137 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
4138
4139 s2 += t1;
4140
4141 /* Renormalize (s1, s2) to (t1, s2) */
4142 t1 = s1 + s2;
4143 s2 = s2 - (t1 - s1);
4144
4145 t2 += s2;
4146
4147 /* Renormalize (t1, t2) */
4148 head_t = t1 + t2;
4149 tail_t = t2 - (head_t - t1);
4150 }
4151 head_sum[1] = head_t;
4152 tail_sum[1] = tail_t;
4153 }
4154 x_elem = tail_x_i[jx];
4155 {
4156 head_prod[0] = (double) a_elem[0] * x_elem;
4157 tail_prod[0] = 0.0;
4158 head_prod[1] = (double) a_elem[1] * x_elem;
4159 tail_prod[1] = 0.0;
4160 }
4161 {
4162 double head_t, tail_t;
4163 double head_a, tail_a;
4164 double head_b, tail_b;
4165 /* Real part */
4166 head_a = head_sum2[0];
4167 tail_a = tail_sum2[0];
4168 head_b = head_prod[0];
4169 tail_b = tail_prod[0];
4170 {
4171 /* Compute double-double = double-double + double-double. */
4172 double bv;
4173 double s1, s2, t1, t2;
4174
4175 /* Add two hi words. */
4176 s1 = head_a + head_b;
4177 bv = s1 - head_a;
4178 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
4179
4180 /* Add two lo words. */
4181 t1 = tail_a + tail_b;
4182 bv = t1 - tail_a;
4183 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
4184
4185 s2 += t1;
4186
4187 /* Renormalize (s1, s2) to (t1, s2) */
4188 t1 = s1 + s2;
4189 s2 = s2 - (t1 - s1);
4190
4191 t2 += s2;
4192
4193 /* Renormalize (t1, t2) */
4194 head_t = t1 + t2;
4195 tail_t = t2 - (head_t - t1);
4196 }
4197 head_sum2[0] = head_t;
4198 tail_sum2[0] = tail_t;
4199 /* Imaginary part */
4200 head_a = head_sum2[1];
4201 tail_a = tail_sum2[1];
4202 head_b = head_prod[1];
4203 tail_b = tail_prod[1];
4204 {
4205 /* Compute double-double = double-double + double-double. */
4206 double bv;
4207 double s1, s2, t1, t2;
4208
4209 /* Add two hi words. */
4210 s1 = head_a + head_b;
4211 bv = s1 - head_a;
4212 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
4213
4214 /* Add two lo words. */
4215 t1 = tail_a + tail_b;
4216 bv = t1 - tail_a;
4217 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
4218
4219 s2 += t1;
4220
4221 /* Renormalize (s1, s2) to (t1, s2) */
4222 t1 = s1 + s2;
4223 s2 = s2 - (t1 - s1);
4224
4225 t2 += s2;
4226
4227 /* Renormalize (t1, t2) */
4228 head_t = t1 + t2;
4229 tail_t = t2 - (head_t - t1);
4230 }
4231 head_sum2[1] = head_t;
4232 tail_sum2[1] = tail_t;
4233 }
4234 aij += incaij;
4235 jx += incx;
4236 }
4237 {
4238 double head_t, tail_t;
4239 double head_a, tail_a;
4240 double head_b, tail_b;
4241 /* Real part */
4242 head_a = head_sum[0];
4243 tail_a = tail_sum[0];
4244 head_b = head_sum2[0];
4245 tail_b = tail_sum2[0];
4246 {
4247 /* Compute double-double = double-double + double-double. */
4248 double bv;
4249 double s1, s2, t1, t2;
4250
4251 /* Add two hi words. */
4252 s1 = head_a + head_b;
4253 bv = s1 - head_a;
4254 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
4255
4256 /* Add two lo words. */
4257 t1 = tail_a + tail_b;
4258 bv = t1 - tail_a;
4259 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
4260
4261 s2 += t1;
4262
4263 /* Renormalize (s1, s2) to (t1, s2) */
4264 t1 = s1 + s2;
4265 s2 = s2 - (t1 - s1);
4266
4267 t2 += s2;
4268
4269 /* Renormalize (t1, t2) */
4270 head_t = t1 + t2;
4271 tail_t = t2 - (head_t - t1);
4272 }
4273 head_sum[0] = head_t;
4274 tail_sum[0] = tail_t;
4275 /* Imaginary part */
4276 head_a = head_sum[1];
4277 tail_a = tail_sum[1];
4278 head_b = head_sum2[1];
4279 tail_b = tail_sum2[1];
4280 {
4281 /* Compute double-double = double-double + double-double. */
4282 double bv;
4283 double s1, s2, t1, t2;
4284
4285 /* Add two hi words. */
4286 s1 = head_a + head_b;
4287 bv = s1 - head_a;
4288 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
4289
4290 /* Add two lo words. */
4291 t1 = tail_a + tail_b;
4292 bv = t1 - tail_a;
4293 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
4294
4295 s2 += t1;
4296
4297 /* Renormalize (s1, s2) to (t1, s2) */
4298 t1 = s1 + s2;
4299 s2 = s2 - (t1 - s1);
4300
4301 t2 += s2;
4302
4303 /* Renormalize (t1, t2) */
4304 head_t = t1 + t2;
4305 tail_t = t2 - (head_t - t1);
4306 }
4307 head_sum[1] = head_t;
4308 tail_sum[1] = tail_t;
4309 }
4310 y_elem[0] = y_i[iy];
4311 y_elem[1] = y_i[iy + 1];
4312 {
4313 double head_e1, tail_e1;
4314 double d1;
4315 double d2;
4316 /* Real part */
4317 d1 = (double) y_elem[0] * beta_i[0];
4318 d2 = (double) -y_elem[1] * beta_i[1];
4319 {
4320 /* Compute double-double = double + double. */
4321 double e, t1, t2;
4322
4323 /* Knuth trick. */
4324 t1 = d1 + d2;
4325 e = t1 - d1;
4326 t2 = ((d2 - e) + (d1 - (t1 - e)));
4327
4328 /* The result is t1 + t2, after normalization. */
4329 head_e1 = t1 + t2;
4330 tail_e1 = t2 - (head_e1 - t1);
4331 }
4332 head_tmp1[0] = head_e1;
4333 tail_tmp1[0] = tail_e1;
4334 /* imaginary part */
4335 d1 = (double) y_elem[0] * beta_i[1];
4336 d2 = (double) y_elem[1] * beta_i[0];
4337 {
4338 /* Compute double-double = double + double. */
4339 double e, t1, t2;
4340
4341 /* Knuth trick. */
4342 t1 = d1 + d2;
4343 e = t1 - d1;
4344 t2 = ((d2 - e) + (d1 - (t1 - e)));
4345
4346 /* The result is t1 + t2, after normalization. */
4347 head_e1 = t1 + t2;
4348 tail_e1 = t2 - (head_e1 - t1);
4349 }
4350 head_tmp1[1] = head_e1;
4351 tail_tmp1[1] = tail_e1;
4352 }
4353 {
4354 double head_t, tail_t;
4355 double head_a, tail_a;
4356 double head_b, tail_b;
4357 /* Real part */
4358 head_a = head_sum[0];
4359 tail_a = tail_sum[0];
4360 head_b = head_tmp1[0];
4361 tail_b = tail_tmp1[0];
4362 {
4363 /* Compute double-double = double-double + double-double. */
4364 double bv;
4365 double s1, s2, t1, t2;
4366
4367 /* Add two hi words. */
4368 s1 = head_a + head_b;
4369 bv = s1 - head_a;
4370 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
4371
4372 /* Add two lo words. */
4373 t1 = tail_a + tail_b;
4374 bv = t1 - tail_a;
4375 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
4376
4377 s2 += t1;
4378
4379 /* Renormalize (s1, s2) to (t1, s2) */
4380 t1 = s1 + s2;
4381 s2 = s2 - (t1 - s1);
4382
4383 t2 += s2;
4384
4385 /* Renormalize (t1, t2) */
4386 head_t = t1 + t2;
4387 tail_t = t2 - (head_t - t1);
4388 }
4389 head_tmp2[0] = head_t;
4390 tail_tmp2[0] = tail_t;
4391 /* Imaginary part */
4392 head_a = head_sum[1];
4393 tail_a = tail_sum[1];
4394 head_b = head_tmp1[1];
4395 tail_b = tail_tmp1[1];
4396 {
4397 /* Compute double-double = double-double + double-double. */
4398 double bv;
4399 double s1, s2, t1, t2;
4400
4401 /* Add two hi words. */
4402 s1 = head_a + head_b;
4403 bv = s1 - head_a;
4404 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
4405
4406 /* Add two lo words. */
4407 t1 = tail_a + tail_b;
4408 bv = t1 - tail_a;
4409 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
4410
4411 s2 += t1;
4412
4413 /* Renormalize (s1, s2) to (t1, s2) */
4414 t1 = s1 + s2;
4415 s2 = s2 - (t1 - s1);
4416
4417 t2 += s2;
4418
4419 /* Renormalize (t1, t2) */
4420 head_t = t1 + t2;
4421 tail_t = t2 - (head_t - t1);
4422 }
4423 head_tmp2[1] = head_t;
4424 tail_tmp2[1] = tail_t;
4425 }
4426 y_i[iy] = head_tmp2[0];
4427 y_i[iy + 1] = head_tmp2[1];
4428 ai += incai;
4429 iy += incy;
4430 }
4431 } else { /* alpha != 1, the most general form:
4432 y = alpha*A*head_x + alpha*A*tail_x + beta*y */
4433 ai = 0;
4434 iy = ky;
4435 for (i = 0; i < leny; i++) {
4436 head_sum[0] = head_sum[1] = tail_sum[0] = tail_sum[1] = 0.0;;
4437 head_sum2[0] = head_sum2[1] = tail_sum2[0] = tail_sum2[1] =
4438 0.0;;
4439 aij = ai;
4440 jx = kx;
4441 for (j = 0; j < lenx; j++) {
4442 a_elem[0] = a_i[aij];
4443 a_elem[1] = a_i[aij + 1];
4444
4445 x_elem = head_x_i[jx];
4446 {
4447 head_prod[0] = (double) a_elem[0] * x_elem;
4448 tail_prod[0] = 0.0;
4449 head_prod[1] = (double) a_elem[1] * x_elem;
4450 tail_prod[1] = 0.0;
4451 }
4452 {
4453 double head_t, tail_t;
4454 double head_a, tail_a;
4455 double head_b, tail_b;
4456 /* Real part */
4457 head_a = head_sum[0];
4458 tail_a = tail_sum[0];
4459 head_b = head_prod[0];
4460 tail_b = tail_prod[0];
4461 {
4462 /* Compute double-double = double-double + double-double. */
4463 double bv;
4464 double s1, s2, t1, t2;
4465
4466 /* Add two hi words. */
4467 s1 = head_a + head_b;
4468 bv = s1 - head_a;
4469 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
4470
4471 /* Add two lo words. */
4472 t1 = tail_a + tail_b;
4473 bv = t1 - tail_a;
4474 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
4475
4476 s2 += t1;
4477
4478 /* Renormalize (s1, s2) to (t1, s2) */
4479 t1 = s1 + s2;
4480 s2 = s2 - (t1 - s1);
4481
4482 t2 += s2;
4483
4484 /* Renormalize (t1, t2) */
4485 head_t = t1 + t2;
4486 tail_t = t2 - (head_t - t1);
4487 }
4488 head_sum[0] = head_t;
4489 tail_sum[0] = tail_t;
4490 /* Imaginary part */
4491 head_a = head_sum[1];
4492 tail_a = tail_sum[1];
4493 head_b = head_prod[1];
4494 tail_b = tail_prod[1];
4495 {
4496 /* Compute double-double = double-double + double-double. */
4497 double bv;
4498 double s1, s2, t1, t2;
4499
4500 /* Add two hi words. */
4501 s1 = head_a + head_b;
4502 bv = s1 - head_a;
4503 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
4504
4505 /* Add two lo words. */
4506 t1 = tail_a + tail_b;
4507 bv = t1 - tail_a;
4508 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
4509
4510 s2 += t1;
4511
4512 /* Renormalize (s1, s2) to (t1, s2) */
4513 t1 = s1 + s2;
4514 s2 = s2 - (t1 - s1);
4515
4516 t2 += s2;
4517
4518 /* Renormalize (t1, t2) */
4519 head_t = t1 + t2;
4520 tail_t = t2 - (head_t - t1);
4521 }
4522 head_sum[1] = head_t;
4523 tail_sum[1] = tail_t;
4524 }
4525 x_elem = tail_x_i[jx];
4526 {
4527 head_prod[0] = (double) a_elem[0] * x_elem;
4528 tail_prod[0] = 0.0;
4529 head_prod[1] = (double) a_elem[1] * x_elem;
4530 tail_prod[1] = 0.0;
4531 }
4532 {
4533 double head_t, tail_t;
4534 double head_a, tail_a;
4535 double head_b, tail_b;
4536 /* Real part */
4537 head_a = head_sum2[0];
4538 tail_a = tail_sum2[0];
4539 head_b = head_prod[0];
4540 tail_b = tail_prod[0];
4541 {
4542 /* Compute double-double = double-double + double-double. */
4543 double bv;
4544 double s1, s2, t1, t2;
4545
4546 /* Add two hi words. */
4547 s1 = head_a + head_b;
4548 bv = s1 - head_a;
4549 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
4550
4551 /* Add two lo words. */
4552 t1 = tail_a + tail_b;
4553 bv = t1 - tail_a;
4554 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
4555
4556 s2 += t1;
4557
4558 /* Renormalize (s1, s2) to (t1, s2) */
4559 t1 = s1 + s2;
4560 s2 = s2 - (t1 - s1);
4561
4562 t2 += s2;
4563
4564 /* Renormalize (t1, t2) */
4565 head_t = t1 + t2;
4566 tail_t = t2 - (head_t - t1);
4567 }
4568 head_sum2[0] = head_t;
4569 tail_sum2[0] = tail_t;
4570 /* Imaginary part */
4571 head_a = head_sum2[1];
4572 tail_a = tail_sum2[1];
4573 head_b = head_prod[1];
4574 tail_b = tail_prod[1];
4575 {
4576 /* Compute double-double = double-double + double-double. */
4577 double bv;
4578 double s1, s2, t1, t2;
4579
4580 /* Add two hi words. */
4581 s1 = head_a + head_b;
4582 bv = s1 - head_a;
4583 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
4584
4585 /* Add two lo words. */
4586 t1 = tail_a + tail_b;
4587 bv = t1 - tail_a;
4588 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
4589
4590 s2 += t1;
4591
4592 /* Renormalize (s1, s2) to (t1, s2) */
4593 t1 = s1 + s2;
4594 s2 = s2 - (t1 - s1);
4595
4596 t2 += s2;
4597
4598 /* Renormalize (t1, t2) */
4599 head_t = t1 + t2;
4600 tail_t = t2 - (head_t - t1);
4601 }
4602 head_sum2[1] = head_t;
4603 tail_sum2[1] = tail_t;
4604 }
4605 aij += incaij;
4606 jx += incx;
4607 }
4608 {
4609 double cd[2];
4610 cd[0] = (double) alpha_i[0];
4611 cd[1] = (double) alpha_i[1];
4612 {
4613 /* Compute complex-extra = complex-extra * complex-double. */
4614 double head_a0, tail_a0;
4615 double head_a1, tail_a1;
4616 double head_t1, tail_t1;
4617 double head_t2, tail_t2;
4618 head_a0 = head_sum[0];
4619 tail_a0 = tail_sum[0];
4620 head_a1 = head_sum[1];
4621 tail_a1 = tail_sum[1];
4622 /* real part */
4623 {
4624 /* Compute double-double = double-double * double. */
4625 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
4626
4627 con = head_a0 * split;
4628 a11 = con - head_a0;
4629 a11 = con - a11;
4630 a21 = head_a0 - a11;
4631 con = cd[0] * split;
4632 b1 = con - cd[0];
4633 b1 = con - b1;
4634 b2 = cd[0] - b1;
4635
4636 c11 = head_a0 * cd[0];
4637 c21 =
4638 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
4639
4640 c2 = tail_a0 * cd[0];
4641 t1 = c11 + c2;
4642 t2 = (c2 - (t1 - c11)) + c21;
4643
4644 head_t1 = t1 + t2;
4645 tail_t1 = t2 - (head_t1 - t1);
4646 }
4647 {
4648 /* Compute double-double = double-double * double. */
4649 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
4650
4651 con = head_a1 * split;
4652 a11 = con - head_a1;
4653 a11 = con - a11;
4654 a21 = head_a1 - a11;
4655 con = cd[1] * split;
4656 b1 = con - cd[1];
4657 b1 = con - b1;
4658 b2 = cd[1] - b1;
4659
4660 c11 = head_a1 * cd[1];
4661 c21 =
4662 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
4663
4664 c2 = tail_a1 * cd[1];
4665 t1 = c11 + c2;
4666 t2 = (c2 - (t1 - c11)) + c21;
4667
4668 head_t2 = t1 + t2;
4669 tail_t2 = t2 - (head_t2 - t1);
4670 }
4671 head_t2 = -head_t2;
4672 tail_t2 = -tail_t2;
4673 {
4674 /* Compute double-double = double-double + double-double. */
4675 double bv;
4676 double s1, s2, t1, t2;
4677
4678 /* Add two hi words. */
4679 s1 = head_t1 + head_t2;
4680 bv = s1 - head_t1;
4681 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
4682
4683 /* Add two lo words. */
4684 t1 = tail_t1 + tail_t2;
4685 bv = t1 - tail_t1;
4686 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
4687
4688 s2 += t1;
4689
4690 /* Renormalize (s1, s2) to (t1, s2) */
4691 t1 = s1 + s2;
4692 s2 = s2 - (t1 - s1);
4693
4694 t2 += s2;
4695
4696 /* Renormalize (t1, t2) */
4697 head_t1 = t1 + t2;
4698 tail_t1 = t2 - (head_t1 - t1);
4699 }
4700 head_tmp1[0] = head_t1;
4701 tail_tmp1[0] = tail_t1;
4702 /* imaginary part */
4703 {
4704 /* Compute double-double = double-double * double. */
4705 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
4706
4707 con = head_a1 * split;
4708 a11 = con - head_a1;
4709 a11 = con - a11;
4710 a21 = head_a1 - a11;
4711 con = cd[0] * split;
4712 b1 = con - cd[0];
4713 b1 = con - b1;
4714 b2 = cd[0] - b1;
4715
4716 c11 = head_a1 * cd[0];
4717 c21 =
4718 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
4719
4720 c2 = tail_a1 * cd[0];
4721 t1 = c11 + c2;
4722 t2 = (c2 - (t1 - c11)) + c21;
4723
4724 head_t1 = t1 + t2;
4725 tail_t1 = t2 - (head_t1 - t1);
4726 }
4727 {
4728 /* Compute double-double = double-double * double. */
4729 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
4730
4731 con = head_a0 * split;
4732 a11 = con - head_a0;
4733 a11 = con - a11;
4734 a21 = head_a0 - a11;
4735 con = cd[1] * split;
4736 b1 = con - cd[1];
4737 b1 = con - b1;
4738 b2 = cd[1] - b1;
4739
4740 c11 = head_a0 * cd[1];
4741 c21 =
4742 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
4743
4744 c2 = tail_a0 * cd[1];
4745 t1 = c11 + c2;
4746 t2 = (c2 - (t1 - c11)) + c21;
4747
4748 head_t2 = t1 + t2;
4749 tail_t2 = t2 - (head_t2 - t1);
4750 }
4751 {
4752 /* Compute double-double = double-double + double-double. */
4753 double bv;
4754 double s1, s2, t1, t2;
4755
4756 /* Add two hi words. */
4757 s1 = head_t1 + head_t2;
4758 bv = s1 - head_t1;
4759 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
4760
4761 /* Add two lo words. */
4762 t1 = tail_t1 + tail_t2;
4763 bv = t1 - tail_t1;
4764 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
4765
4766 s2 += t1;
4767
4768 /* Renormalize (s1, s2) to (t1, s2) */
4769 t1 = s1 + s2;
4770 s2 = s2 - (t1 - s1);
4771
4772 t2 += s2;
4773
4774 /* Renormalize (t1, t2) */
4775 head_t1 = t1 + t2;
4776 tail_t1 = t2 - (head_t1 - t1);
4777 }
4778 head_tmp1[1] = head_t1;
4779 tail_tmp1[1] = tail_t1;
4780 }
4781
4782 }
4783 {
4784 double cd[2];
4785 cd[0] = (double) alpha_i[0];
4786 cd[1] = (double) alpha_i[1];
4787 {
4788 /* Compute complex-extra = complex-extra * complex-double. */
4789 double head_a0, tail_a0;
4790 double head_a1, tail_a1;
4791 double head_t1, tail_t1;
4792 double head_t2, tail_t2;
4793 head_a0 = head_sum2[0];
4794 tail_a0 = tail_sum2[0];
4795 head_a1 = head_sum2[1];
4796 tail_a1 = tail_sum2[1];
4797 /* real part */
4798 {
4799 /* Compute double-double = double-double * double. */
4800 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
4801
4802 con = head_a0 * split;
4803 a11 = con - head_a0;
4804 a11 = con - a11;
4805 a21 = head_a0 - a11;
4806 con = cd[0] * split;
4807 b1 = con - cd[0];
4808 b1 = con - b1;
4809 b2 = cd[0] - b1;
4810
4811 c11 = head_a0 * cd[0];
4812 c21 =
4813 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
4814
4815 c2 = tail_a0 * cd[0];
4816 t1 = c11 + c2;
4817 t2 = (c2 - (t1 - c11)) + c21;
4818
4819 head_t1 = t1 + t2;
4820 tail_t1 = t2 - (head_t1 - t1);
4821 }
4822 {
4823 /* Compute double-double = double-double * double. */
4824 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
4825
4826 con = head_a1 * split;
4827 a11 = con - head_a1;
4828 a11 = con - a11;
4829 a21 = head_a1 - a11;
4830 con = cd[1] * split;
4831 b1 = con - cd[1];
4832 b1 = con - b1;
4833 b2 = cd[1] - b1;
4834
4835 c11 = head_a1 * cd[1];
4836 c21 =
4837 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
4838
4839 c2 = tail_a1 * cd[1];
4840 t1 = c11 + c2;
4841 t2 = (c2 - (t1 - c11)) + c21;
4842
4843 head_t2 = t1 + t2;
4844 tail_t2 = t2 - (head_t2 - t1);
4845 }
4846 head_t2 = -head_t2;
4847 tail_t2 = -tail_t2;
4848 {
4849 /* Compute double-double = double-double + double-double. */
4850 double bv;
4851 double s1, s2, t1, t2;
4852
4853 /* Add two hi words. */
4854 s1 = head_t1 + head_t2;
4855 bv = s1 - head_t1;
4856 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
4857
4858 /* Add two lo words. */
4859 t1 = tail_t1 + tail_t2;
4860 bv = t1 - tail_t1;
4861 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
4862
4863 s2 += t1;
4864
4865 /* Renormalize (s1, s2) to (t1, s2) */
4866 t1 = s1 + s2;
4867 s2 = s2 - (t1 - s1);
4868
4869 t2 += s2;
4870
4871 /* Renormalize (t1, t2) */
4872 head_t1 = t1 + t2;
4873 tail_t1 = t2 - (head_t1 - t1);
4874 }
4875 head_tmp2[0] = head_t1;
4876 tail_tmp2[0] = tail_t1;
4877 /* imaginary part */
4878 {
4879 /* Compute double-double = double-double * double. */
4880 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
4881
4882 con = head_a1 * split;
4883 a11 = con - head_a1;
4884 a11 = con - a11;
4885 a21 = head_a1 - a11;
4886 con = cd[0] * split;
4887 b1 = con - cd[0];
4888 b1 = con - b1;
4889 b2 = cd[0] - b1;
4890
4891 c11 = head_a1 * cd[0];
4892 c21 =
4893 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
4894
4895 c2 = tail_a1 * cd[0];
4896 t1 = c11 + c2;
4897 t2 = (c2 - (t1 - c11)) + c21;
4898
4899 head_t1 = t1 + t2;
4900 tail_t1 = t2 - (head_t1 - t1);
4901 }
4902 {
4903 /* Compute double-double = double-double * double. */
4904 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
4905
4906 con = head_a0 * split;
4907 a11 = con - head_a0;
4908 a11 = con - a11;
4909 a21 = head_a0 - a11;
4910 con = cd[1] * split;
4911 b1 = con - cd[1];
4912 b1 = con - b1;
4913 b2 = cd[1] - b1;
4914
4915 c11 = head_a0 * cd[1];
4916 c21 =
4917 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
4918
4919 c2 = tail_a0 * cd[1];
4920 t1 = c11 + c2;
4921 t2 = (c2 - (t1 - c11)) + c21;
4922
4923 head_t2 = t1 + t2;
4924 tail_t2 = t2 - (head_t2 - t1);
4925 }
4926 {
4927 /* Compute double-double = double-double + double-double. */
4928 double bv;
4929 double s1, s2, t1, t2;
4930
4931 /* Add two hi words. */
4932 s1 = head_t1 + head_t2;
4933 bv = s1 - head_t1;
4934 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
4935
4936 /* Add two lo words. */
4937 t1 = tail_t1 + tail_t2;
4938 bv = t1 - tail_t1;
4939 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
4940
4941 s2 += t1;
4942
4943 /* Renormalize (s1, s2) to (t1, s2) */
4944 t1 = s1 + s2;
4945 s2 = s2 - (t1 - s1);
4946
4947 t2 += s2;
4948
4949 /* Renormalize (t1, t2) */
4950 head_t1 = t1 + t2;
4951 tail_t1 = t2 - (head_t1 - t1);
4952 }
4953 head_tmp2[1] = head_t1;
4954 tail_tmp2[1] = tail_t1;
4955 }
4956
4957 }
4958 {
4959 double head_t, tail_t;
4960 double head_a, tail_a;
4961 double head_b, tail_b;
4962 /* Real part */
4963 head_a = head_tmp1[0];
4964 tail_a = tail_tmp1[0];
4965 head_b = head_tmp2[0];
4966 tail_b = tail_tmp2[0];
4967 {
4968 /* Compute double-double = double-double + double-double. */
4969 double bv;
4970 double s1, s2, t1, t2;
4971
4972 /* Add two hi words. */
4973 s1 = head_a + head_b;
4974 bv = s1 - head_a;
4975 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
4976
4977 /* Add two lo words. */
4978 t1 = tail_a + tail_b;
4979 bv = t1 - tail_a;
4980 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
4981
4982 s2 += t1;
4983
4984 /* Renormalize (s1, s2) to (t1, s2) */
4985 t1 = s1 + s2;
4986 s2 = s2 - (t1 - s1);
4987
4988 t2 += s2;
4989
4990 /* Renormalize (t1, t2) */
4991 head_t = t1 + t2;
4992 tail_t = t2 - (head_t - t1);
4993 }
4994 head_tmp1[0] = head_t;
4995 tail_tmp1[0] = tail_t;
4996 /* Imaginary part */
4997 head_a = head_tmp1[1];
4998 tail_a = tail_tmp1[1];
4999 head_b = head_tmp2[1];
5000 tail_b = tail_tmp2[1];
5001 {
5002 /* Compute double-double = double-double + double-double. */
5003 double bv;
5004 double s1, s2, t1, t2;
5005
5006 /* Add two hi words. */
5007 s1 = head_a + head_b;
5008 bv = s1 - head_a;
5009 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
5010
5011 /* Add two lo words. */
5012 t1 = tail_a + tail_b;
5013 bv = t1 - tail_a;
5014 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
5015
5016 s2 += t1;
5017
5018 /* Renormalize (s1, s2) to (t1, s2) */
5019 t1 = s1 + s2;
5020 s2 = s2 - (t1 - s1);
5021
5022 t2 += s2;
5023
5024 /* Renormalize (t1, t2) */
5025 head_t = t1 + t2;
5026 tail_t = t2 - (head_t - t1);
5027 }
5028 head_tmp1[1] = head_t;
5029 tail_tmp1[1] = tail_t;
5030 }
5031 y_elem[0] = y_i[iy];
5032 y_elem[1] = y_i[iy + 1];
5033 {
5034 double head_e1, tail_e1;
5035 double d1;
5036 double d2;
5037 /* Real part */
5038 d1 = (double) y_elem[0] * beta_i[0];
5039 d2 = (double) -y_elem[1] * beta_i[1];
5040 {
5041 /* Compute double-double = double + double. */
5042 double e, t1, t2;
5043
5044 /* Knuth trick. */
5045 t1 = d1 + d2;
5046 e = t1 - d1;
5047 t2 = ((d2 - e) + (d1 - (t1 - e)));
5048
5049 /* The result is t1 + t2, after normalization. */
5050 head_e1 = t1 + t2;
5051 tail_e1 = t2 - (head_e1 - t1);
5052 }
5053 head_tmp2[0] = head_e1;
5054 tail_tmp2[0] = tail_e1;
5055 /* imaginary part */
5056 d1 = (double) y_elem[0] * beta_i[1];
5057 d2 = (double) y_elem[1] * beta_i[0];
5058 {
5059 /* Compute double-double = double + double. */
5060 double e, t1, t2;
5061
5062 /* Knuth trick. */
5063 t1 = d1 + d2;
5064 e = t1 - d1;
5065 t2 = ((d2 - e) + (d1 - (t1 - e)));
5066
5067 /* The result is t1 + t2, after normalization. */
5068 head_e1 = t1 + t2;
5069 tail_e1 = t2 - (head_e1 - t1);
5070 }
5071 head_tmp2[1] = head_e1;
5072 tail_tmp2[1] = tail_e1;
5073 }
5074 {
5075 double head_t, tail_t;
5076 double head_a, tail_a;
5077 double head_b, tail_b;
5078 /* Real part */
5079 head_a = head_tmp1[0];
5080 tail_a = tail_tmp1[0];
5081 head_b = head_tmp2[0];
5082 tail_b = tail_tmp2[0];
5083 {
5084 /* Compute double-double = double-double + double-double. */
5085 double bv;
5086 double s1, s2, t1, t2;
5087
5088 /* Add two hi words. */
5089 s1 = head_a + head_b;
5090 bv = s1 - head_a;
5091 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
5092
5093 /* Add two lo words. */
5094 t1 = tail_a + tail_b;
5095 bv = t1 - tail_a;
5096 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
5097
5098 s2 += t1;
5099
5100 /* Renormalize (s1, s2) to (t1, s2) */
5101 t1 = s1 + s2;
5102 s2 = s2 - (t1 - s1);
5103
5104 t2 += s2;
5105
5106 /* Renormalize (t1, t2) */
5107 head_t = t1 + t2;
5108 tail_t = t2 - (head_t - t1);
5109 }
5110 head_tmp1[0] = head_t;
5111 tail_tmp1[0] = tail_t;
5112 /* Imaginary part */
5113 head_a = head_tmp1[1];
5114 tail_a = tail_tmp1[1];
5115 head_b = head_tmp2[1];
5116 tail_b = tail_tmp2[1];
5117 {
5118 /* Compute double-double = double-double + double-double. */
5119 double bv;
5120 double s1, s2, t1, t2;
5121
5122 /* Add two hi words. */
5123 s1 = head_a + head_b;
5124 bv = s1 - head_a;
5125 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
5126
5127 /* Add two lo words. */
5128 t1 = tail_a + tail_b;
5129 bv = t1 - tail_a;
5130 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
5131
5132 s2 += t1;
5133
5134 /* Renormalize (s1, s2) to (t1, s2) */
5135 t1 = s1 + s2;
5136 s2 = s2 - (t1 - s1);
5137
5138 t2 += s2;
5139
5140 /* Renormalize (t1, t2) */
5141 head_t = t1 + t2;
5142 tail_t = t2 - (head_t - t1);
5143 }
5144 head_tmp1[1] = head_t;
5145 tail_tmp1[1] = tail_t;
5146 }
5147 y_i[iy] = head_tmp1[0];
5148 y_i[iy + 1] = head_tmp1[1];
5149 ai += incai;
5150 iy += incy;
5151 }
5152 }
5153 }
5154
5155 }
5156 }
5157
5158 FPU_FIX_STOP;
5159 }
5160 break;
5161 }
5162 }
5163