1 #include "blas_extended.h"
2 #include "blas_extended_private.h"
BLAS_cgemv_x(enum blas_order_type order,enum blas_trans_type trans,int m,int n,const void * alpha,const void * a,int lda,const void * x,int incx,const void * beta,void * y,int incy,enum blas_prec_type prec)3 void BLAS_cgemv_x(enum blas_order_type order, enum blas_trans_type trans,
4 int m, int n, const void *alpha, const void *a, int lda,
5 const void *x, int incx, const void *beta, void *y,
6 int incy, enum blas_prec_type prec)
7
8 /*
9 * Purpose
10 * =======
11 *
12 * Computes y = alpha * A * x + beta * y, where A is a general matrix.
13 *
14 * Arguments
15 * =========
16 *
17 * order (input) blas_order_type
18 * Order of AP; row or column major
19 *
20 * trans (input) blas_trans_type
21 * Transpose of AB; no trans,
22 * trans, or conjugate trans
23 *
24 * m (input) int
25 * Dimension of AB
26 *
27 * n (input) int
28 * Dimension of AB and the length of vector x
29 *
30 * alpha (input) const void*
31 *
32 * A (input) const void*
33 *
34 * lda (input) int
35 * Leading dimension of A
36 *
37 * x (input) const void*
38 *
39 * incx (input) int
40 * The stride for vector x.
41 *
42 * beta (input) const void*
43 *
44 * y (input/output) void*
45 *
46 * incy (input) int
47 * The stride for vector y.
48 *
49 * prec (input) enum blas_prec_type
50 * Specifies the internal precision to be used.
51 * = blas_prec_single: single precision.
52 * = blas_prec_double: double precision.
53 * = blas_prec_extra : anything at least 1.5 times as accurate
54 * than double, and wider than 80-bits.
55 * We use double-double in our implementation.
56 *
57 */
58 {
59 static const char routine_name[] = "BLAS_cgemv_x";
60 switch (prec) {
61 case blas_prec_single:{
62
63 int i, j;
64 int iy, jx, kx, ky;
65 int lenx, leny;
66 int ai, aij;
67 int incai, incaij;
68
69 const float *a_i = (float *) a;
70 const float *x_i = (float *) x;
71 float *y_i = (float *) y;
72 float *alpha_i = (float *) alpha;
73 float *beta_i = (float *) beta;
74 float a_elem[2];
75 float x_elem[2];
76 float y_elem[2];
77 float prod[2];
78 float sum[2];
79 float tmp1[2];
80 float tmp2[2];
81
82
83 /* all error calls */
84 if (m < 0)
85 BLAS_error(routine_name, -3, m, 0);
86 else if (n <= 0)
87 BLAS_error(routine_name, -4, n, 0);
88 else if (incx == 0)
89 BLAS_error(routine_name, -9, incx, 0);
90 else if (incy == 0)
91 BLAS_error(routine_name, -12, incy, 0);
92
93 if ((order == blas_rowmajor) && (trans == blas_no_trans)) {
94 lenx = n;
95 leny = m;
96 incai = lda;
97 incaij = 1;
98 } else if ((order == blas_rowmajor) && (trans != blas_no_trans)) {
99 lenx = m;
100 leny = n;
101 incai = 1;
102 incaij = lda;
103 } else if ((order == blas_colmajor) && (trans == blas_no_trans)) {
104 lenx = n;
105 leny = m;
106 incai = 1;
107 incaij = lda;
108 } else { /* colmajor and blas_trans */
109 lenx = m;
110 leny = n;
111 incai = lda;
112 incaij = 1;
113 }
114 if ((order == blas_colmajor && lda < m) ||
115 (order == blas_rowmajor && lda < n))
116 BLAS_error(routine_name, -7, lda, NULL);
117
118
119
120 incx *= 2;
121 incy *= 2;
122 incai *= 2;
123 incaij *= 2;
124
125 if (incx > 0)
126 kx = 0;
127 else
128 kx = (1 - lenx) * incx;
129 if (incy > 0)
130 ky = 0;
131 else
132 ky = (1 - leny) * incy;
133
134 /* No extra-precision needed for alpha = 0 */
135 if (alpha_i[0] == 0.0 && alpha_i[1] == 0.0) {
136 if (beta_i[0] == 0.0 && beta_i[1] == 0.0) {
137 iy = ky;
138 for (i = 0; i < leny; i++) {
139 y_i[iy] = 0.0;
140 y_i[iy + 1] = 0.0;
141 iy += incy;
142 }
143 } else if (!(beta_i[0] == 0.0 && beta_i[1] == 0.0)) {
144 iy = ky;
145 for (i = 0; i < leny; i++) {
146 y_elem[0] = y_i[iy];
147 y_elem[1] = y_i[iy + 1];
148 {
149 tmp1[0] = y_elem[0] * beta_i[0] - y_elem[1] * beta_i[1];
150 tmp1[1] = y_elem[0] * beta_i[1] + y_elem[1] * beta_i[0];
151 }
152
153 y_i[iy] = tmp1[0];
154 y_i[iy + 1] = tmp1[1];
155 iy += incy;
156 }
157 }
158 } else {
159 if (trans == blas_conj_trans) {
160
161 /* if beta = 0, we can save m multiplies: y = alpha*A*x */
162 if (beta_i[0] == 0.0 && beta_i[1] == 0.0) {
163 /* save m more multiplies if alpha = 1 */
164 if ((alpha_i[0] == 1.0 && alpha_i[1] == 0.0)) {
165 ai = 0;
166 iy = ky;
167 for (i = 0; i < leny; i++) {
168 sum[0] = sum[1] = 0.0;
169 aij = ai;
170 jx = kx;
171 for (j = 0; j < lenx; j++) {
172 a_elem[0] = a_i[aij];
173 a_elem[1] = a_i[aij + 1];
174 a_elem[1] = -a_elem[1];
175 x_elem[0] = x_i[jx];
176 x_elem[1] = x_i[jx + 1];
177 {
178 prod[0] = a_elem[0] * x_elem[0] - a_elem[1] * x_elem[1];
179 prod[1] = a_elem[0] * x_elem[1] + a_elem[1] * x_elem[0];
180 }
181
182 sum[0] = sum[0] + prod[0];
183 sum[1] = sum[1] + prod[1];
184 aij += incaij;
185 jx += incx;
186 }
187 y_i[iy] = sum[0];
188 y_i[iy + 1] = sum[1];
189 ai += incai;
190 iy += incy;
191 }
192 } else {
193 ai = 0;
194 iy = ky;
195 for (i = 0; i < leny; i++) {
196 sum[0] = sum[1] = 0.0;
197 aij = ai;
198 jx = kx;
199 for (j = 0; j < lenx; j++) {
200 a_elem[0] = a_i[aij];
201 a_elem[1] = a_i[aij + 1];
202 a_elem[1] = -a_elem[1];
203 x_elem[0] = x_i[jx];
204 x_elem[1] = x_i[jx + 1];
205 {
206 prod[0] = a_elem[0] * x_elem[0] - a_elem[1] * x_elem[1];
207 prod[1] = a_elem[0] * x_elem[1] + a_elem[1] * x_elem[0];
208 }
209
210 sum[0] = sum[0] + prod[0];
211 sum[1] = sum[1] + prod[1];
212 aij += incaij;
213 jx += incx;
214 }
215 {
216 tmp1[0] = sum[0] * alpha_i[0] - sum[1] * alpha_i[1];
217 tmp1[1] = sum[0] * alpha_i[1] + sum[1] * alpha_i[0];
218 }
219
220 y_i[iy] = tmp1[0];
221 y_i[iy + 1] = tmp1[1];
222 ai += incai;
223 iy += incy;
224 }
225 }
226 } else {
227 /* the most general form, y = alpha*A*x + beta*y */
228 ai = 0;
229 iy = ky;
230 for (i = 0; i < leny; i++) {
231 sum[0] = sum[1] = 0.0;;
232 aij = ai;
233 jx = kx;
234 for (j = 0; j < lenx; j++) {
235 a_elem[0] = a_i[aij];
236 a_elem[1] = a_i[aij + 1];
237 a_elem[1] = -a_elem[1];
238 x_elem[0] = x_i[jx];
239 x_elem[1] = x_i[jx + 1];
240 {
241 prod[0] = a_elem[0] * x_elem[0] - a_elem[1] * x_elem[1];
242 prod[1] = a_elem[0] * x_elem[1] + a_elem[1] * x_elem[0];
243 }
244
245 sum[0] = sum[0] + prod[0];
246 sum[1] = sum[1] + prod[1];
247 aij += incaij;
248 jx += incx;
249 }
250 {
251 tmp1[0] = sum[0] * alpha_i[0] - sum[1] * alpha_i[1];
252 tmp1[1] = sum[0] * alpha_i[1] + sum[1] * alpha_i[0];
253 }
254
255 y_elem[0] = y_i[iy];
256 y_elem[1] = y_i[iy + 1];
257 {
258 tmp2[0] = y_elem[0] * beta_i[0] - y_elem[1] * beta_i[1];
259 tmp2[1] = y_elem[0] * beta_i[1] + y_elem[1] * beta_i[0];
260 }
261
262 tmp1[0] = tmp1[0] + tmp2[0];
263 tmp1[1] = tmp1[1] + tmp2[1];
264 y_i[iy] = tmp1[0];
265 y_i[iy + 1] = tmp1[1];
266 ai += incai;
267 iy += incy;
268 }
269 }
270
271 } else {
272
273 /* if beta = 0, we can save m multiplies: y = alpha*A*x */
274 if (beta_i[0] == 0.0 && beta_i[1] == 0.0) {
275 /* save m more multiplies if alpha = 1 */
276 if ((alpha_i[0] == 1.0 && alpha_i[1] == 0.0)) {
277 ai = 0;
278 iy = ky;
279 for (i = 0; i < leny; i++) {
280 sum[0] = sum[1] = 0.0;
281 aij = ai;
282 jx = kx;
283 for (j = 0; j < lenx; j++) {
284 a_elem[0] = a_i[aij];
285 a_elem[1] = a_i[aij + 1];
286
287 x_elem[0] = x_i[jx];
288 x_elem[1] = x_i[jx + 1];
289 {
290 prod[0] = a_elem[0] * x_elem[0] - a_elem[1] * x_elem[1];
291 prod[1] = a_elem[0] * x_elem[1] + a_elem[1] * x_elem[0];
292 }
293
294 sum[0] = sum[0] + prod[0];
295 sum[1] = sum[1] + prod[1];
296 aij += incaij;
297 jx += incx;
298 }
299 y_i[iy] = sum[0];
300 y_i[iy + 1] = sum[1];
301 ai += incai;
302 iy += incy;
303 }
304 } else {
305 ai = 0;
306 iy = ky;
307 for (i = 0; i < leny; i++) {
308 sum[0] = sum[1] = 0.0;
309 aij = ai;
310 jx = kx;
311 for (j = 0; j < lenx; j++) {
312 a_elem[0] = a_i[aij];
313 a_elem[1] = a_i[aij + 1];
314
315 x_elem[0] = x_i[jx];
316 x_elem[1] = x_i[jx + 1];
317 {
318 prod[0] = a_elem[0] * x_elem[0] - a_elem[1] * x_elem[1];
319 prod[1] = a_elem[0] * x_elem[1] + a_elem[1] * x_elem[0];
320 }
321
322 sum[0] = sum[0] + prod[0];
323 sum[1] = sum[1] + prod[1];
324 aij += incaij;
325 jx += incx;
326 }
327 {
328 tmp1[0] = sum[0] * alpha_i[0] - sum[1] * alpha_i[1];
329 tmp1[1] = sum[0] * alpha_i[1] + sum[1] * alpha_i[0];
330 }
331
332 y_i[iy] = tmp1[0];
333 y_i[iy + 1] = tmp1[1];
334 ai += incai;
335 iy += incy;
336 }
337 }
338 } else {
339 /* the most general form, y = alpha*A*x + beta*y */
340 ai = 0;
341 iy = ky;
342 for (i = 0; i < leny; i++) {
343 sum[0] = sum[1] = 0.0;;
344 aij = ai;
345 jx = kx;
346 for (j = 0; j < lenx; j++) {
347 a_elem[0] = a_i[aij];
348 a_elem[1] = a_i[aij + 1];
349
350 x_elem[0] = x_i[jx];
351 x_elem[1] = x_i[jx + 1];
352 {
353 prod[0] = a_elem[0] * x_elem[0] - a_elem[1] * x_elem[1];
354 prod[1] = a_elem[0] * x_elem[1] + a_elem[1] * x_elem[0];
355 }
356
357 sum[0] = sum[0] + prod[0];
358 sum[1] = sum[1] + prod[1];
359 aij += incaij;
360 jx += incx;
361 }
362 {
363 tmp1[0] = sum[0] * alpha_i[0] - sum[1] * alpha_i[1];
364 tmp1[1] = sum[0] * alpha_i[1] + sum[1] * alpha_i[0];
365 }
366
367 y_elem[0] = y_i[iy];
368 y_elem[1] = y_i[iy + 1];
369 {
370 tmp2[0] = y_elem[0] * beta_i[0] - y_elem[1] * beta_i[1];
371 tmp2[1] = y_elem[0] * beta_i[1] + y_elem[1] * beta_i[0];
372 }
373
374 tmp1[0] = tmp1[0] + tmp2[0];
375 tmp1[1] = tmp1[1] + tmp2[1];
376 y_i[iy] = tmp1[0];
377 y_i[iy + 1] = tmp1[1];
378 ai += incai;
379 iy += incy;
380 }
381 }
382
383 }
384 }
385
386
387
388 break;
389 }
390 case blas_prec_double:
391 case blas_prec_indigenous:{
392
393 int i, j;
394 int iy, jx, kx, ky;
395 int lenx, leny;
396 int ai, aij;
397 int incai, incaij;
398
399 const float *a_i = (float *) a;
400 const float *x_i = (float *) x;
401 float *y_i = (float *) y;
402 float *alpha_i = (float *) alpha;
403 float *beta_i = (float *) beta;
404 float a_elem[2];
405 float x_elem[2];
406 float y_elem[2];
407 double prod[2];
408 double sum[2];
409 double tmp1[2];
410 double tmp2[2];
411
412
413 /* all error calls */
414 if (m < 0)
415 BLAS_error(routine_name, -3, m, 0);
416 else if (n <= 0)
417 BLAS_error(routine_name, -4, n, 0);
418 else if (incx == 0)
419 BLAS_error(routine_name, -9, incx, 0);
420 else if (incy == 0)
421 BLAS_error(routine_name, -12, incy, 0);
422
423 if ((order == blas_rowmajor) && (trans == blas_no_trans)) {
424 lenx = n;
425 leny = m;
426 incai = lda;
427 incaij = 1;
428 } else if ((order == blas_rowmajor) && (trans != blas_no_trans)) {
429 lenx = m;
430 leny = n;
431 incai = 1;
432 incaij = lda;
433 } else if ((order == blas_colmajor) && (trans == blas_no_trans)) {
434 lenx = n;
435 leny = m;
436 incai = 1;
437 incaij = lda;
438 } else { /* colmajor and blas_trans */
439 lenx = m;
440 leny = n;
441 incai = lda;
442 incaij = 1;
443 }
444 if ((order == blas_colmajor && lda < m) ||
445 (order == blas_rowmajor && lda < n))
446 BLAS_error(routine_name, -7, lda, NULL);
447
448
449
450 incx *= 2;
451 incy *= 2;
452 incai *= 2;
453 incaij *= 2;
454
455 if (incx > 0)
456 kx = 0;
457 else
458 kx = (1 - lenx) * incx;
459 if (incy > 0)
460 ky = 0;
461 else
462 ky = (1 - leny) * incy;
463
464 /* No extra-precision needed for alpha = 0 */
465 if (alpha_i[0] == 0.0 && alpha_i[1] == 0.0) {
466 if (beta_i[0] == 0.0 && beta_i[1] == 0.0) {
467 iy = ky;
468 for (i = 0; i < leny; i++) {
469 y_i[iy] = 0.0;
470 y_i[iy + 1] = 0.0;
471 iy += incy;
472 }
473 } else if (!(beta_i[0] == 0.0 && beta_i[1] == 0.0)) {
474 iy = ky;
475 for (i = 0; i < leny; i++) {
476 y_elem[0] = y_i[iy];
477 y_elem[1] = y_i[iy + 1];
478 {
479 tmp1[0] =
480 (double) y_elem[0] * beta_i[0] -
481 (double) y_elem[1] * beta_i[1];
482 tmp1[1] =
483 (double) y_elem[0] * beta_i[1] +
484 (double) y_elem[1] * beta_i[0];
485 }
486 y_i[iy] = tmp1[0];
487 y_i[iy + 1] = tmp1[1];
488 iy += incy;
489 }
490 }
491 } else {
492 if (trans == blas_conj_trans) {
493
494 /* if beta = 0, we can save m multiplies: y = alpha*A*x */
495 if (beta_i[0] == 0.0 && beta_i[1] == 0.0) {
496 /* save m more multiplies if alpha = 1 */
497 if ((alpha_i[0] == 1.0 && alpha_i[1] == 0.0)) {
498 ai = 0;
499 iy = ky;
500 for (i = 0; i < leny; i++) {
501 sum[0] = sum[1] = 0.0;
502 aij = ai;
503 jx = kx;
504 for (j = 0; j < lenx; j++) {
505 a_elem[0] = a_i[aij];
506 a_elem[1] = a_i[aij + 1];
507 a_elem[1] = -a_elem[1];
508 x_elem[0] = x_i[jx];
509 x_elem[1] = x_i[jx + 1];
510 {
511 prod[0] =
512 (double) a_elem[0] * x_elem[0] -
513 (double) a_elem[1] * x_elem[1];
514 prod[1] =
515 (double) a_elem[0] * x_elem[1] +
516 (double) a_elem[1] * x_elem[0];
517 }
518 sum[0] = sum[0] + prod[0];
519 sum[1] = sum[1] + prod[1];
520 aij += incaij;
521 jx += incx;
522 }
523 y_i[iy] = sum[0];
524 y_i[iy + 1] = sum[1];
525 ai += incai;
526 iy += incy;
527 }
528 } else {
529 ai = 0;
530 iy = ky;
531 for (i = 0; i < leny; i++) {
532 sum[0] = sum[1] = 0.0;
533 aij = ai;
534 jx = kx;
535 for (j = 0; j < lenx; j++) {
536 a_elem[0] = a_i[aij];
537 a_elem[1] = a_i[aij + 1];
538 a_elem[1] = -a_elem[1];
539 x_elem[0] = x_i[jx];
540 x_elem[1] = x_i[jx + 1];
541 {
542 prod[0] =
543 (double) a_elem[0] * x_elem[0] -
544 (double) a_elem[1] * x_elem[1];
545 prod[1] =
546 (double) a_elem[0] * x_elem[1] +
547 (double) a_elem[1] * x_elem[0];
548 }
549 sum[0] = sum[0] + prod[0];
550 sum[1] = sum[1] + prod[1];
551 aij += incaij;
552 jx += incx;
553 }
554 {
555 tmp1[0] =
556 (double) sum[0] * alpha_i[0] -
557 (double) sum[1] * alpha_i[1];
558 tmp1[1] =
559 (double) sum[0] * alpha_i[1] +
560 (double) sum[1] * alpha_i[0];
561 }
562 y_i[iy] = tmp1[0];
563 y_i[iy + 1] = tmp1[1];
564 ai += incai;
565 iy += incy;
566 }
567 }
568 } else {
569 /* the most general form, y = alpha*A*x + beta*y */
570 ai = 0;
571 iy = ky;
572 for (i = 0; i < leny; i++) {
573 sum[0] = sum[1] = 0.0;;
574 aij = ai;
575 jx = kx;
576 for (j = 0; j < lenx; j++) {
577 a_elem[0] = a_i[aij];
578 a_elem[1] = a_i[aij + 1];
579 a_elem[1] = -a_elem[1];
580 x_elem[0] = x_i[jx];
581 x_elem[1] = x_i[jx + 1];
582 {
583 prod[0] =
584 (double) a_elem[0] * x_elem[0] -
585 (double) a_elem[1] * x_elem[1];
586 prod[1] =
587 (double) a_elem[0] * x_elem[1] +
588 (double) a_elem[1] * x_elem[0];
589 }
590 sum[0] = sum[0] + prod[0];
591 sum[1] = sum[1] + prod[1];
592 aij += incaij;
593 jx += incx;
594 }
595 {
596 tmp1[0] =
597 (double) sum[0] * alpha_i[0] - (double) sum[1] * alpha_i[1];
598 tmp1[1] =
599 (double) sum[0] * alpha_i[1] + (double) sum[1] * alpha_i[0];
600 }
601 y_elem[0] = y_i[iy];
602 y_elem[1] = y_i[iy + 1];
603 {
604 tmp2[0] =
605 (double) y_elem[0] * beta_i[0] -
606 (double) y_elem[1] * beta_i[1];
607 tmp2[1] =
608 (double) y_elem[0] * beta_i[1] +
609 (double) y_elem[1] * beta_i[0];
610 }
611 tmp1[0] = tmp1[0] + tmp2[0];
612 tmp1[1] = tmp1[1] + tmp2[1];
613 y_i[iy] = tmp1[0];
614 y_i[iy + 1] = tmp1[1];
615 ai += incai;
616 iy += incy;
617 }
618 }
619
620 } else {
621
622 /* if beta = 0, we can save m multiplies: y = alpha*A*x */
623 if (beta_i[0] == 0.0 && beta_i[1] == 0.0) {
624 /* save m more multiplies if alpha = 1 */
625 if ((alpha_i[0] == 1.0 && alpha_i[1] == 0.0)) {
626 ai = 0;
627 iy = ky;
628 for (i = 0; i < leny; i++) {
629 sum[0] = sum[1] = 0.0;
630 aij = ai;
631 jx = kx;
632 for (j = 0; j < lenx; j++) {
633 a_elem[0] = a_i[aij];
634 a_elem[1] = a_i[aij + 1];
635
636 x_elem[0] = x_i[jx];
637 x_elem[1] = x_i[jx + 1];
638 {
639 prod[0] =
640 (double) a_elem[0] * x_elem[0] -
641 (double) a_elem[1] * x_elem[1];
642 prod[1] =
643 (double) a_elem[0] * x_elem[1] +
644 (double) a_elem[1] * x_elem[0];
645 }
646 sum[0] = sum[0] + prod[0];
647 sum[1] = sum[1] + prod[1];
648 aij += incaij;
649 jx += incx;
650 }
651 y_i[iy] = sum[0];
652 y_i[iy + 1] = sum[1];
653 ai += incai;
654 iy += incy;
655 }
656 } else {
657 ai = 0;
658 iy = ky;
659 for (i = 0; i < leny; i++) {
660 sum[0] = sum[1] = 0.0;
661 aij = ai;
662 jx = kx;
663 for (j = 0; j < lenx; j++) {
664 a_elem[0] = a_i[aij];
665 a_elem[1] = a_i[aij + 1];
666
667 x_elem[0] = x_i[jx];
668 x_elem[1] = x_i[jx + 1];
669 {
670 prod[0] =
671 (double) a_elem[0] * x_elem[0] -
672 (double) a_elem[1] * x_elem[1];
673 prod[1] =
674 (double) a_elem[0] * x_elem[1] +
675 (double) a_elem[1] * x_elem[0];
676 }
677 sum[0] = sum[0] + prod[0];
678 sum[1] = sum[1] + prod[1];
679 aij += incaij;
680 jx += incx;
681 }
682 {
683 tmp1[0] =
684 (double) sum[0] * alpha_i[0] -
685 (double) sum[1] * alpha_i[1];
686 tmp1[1] =
687 (double) sum[0] * alpha_i[1] +
688 (double) sum[1] * alpha_i[0];
689 }
690 y_i[iy] = tmp1[0];
691 y_i[iy + 1] = tmp1[1];
692 ai += incai;
693 iy += incy;
694 }
695 }
696 } else {
697 /* the most general form, y = alpha*A*x + beta*y */
698 ai = 0;
699 iy = ky;
700 for (i = 0; i < leny; i++) {
701 sum[0] = sum[1] = 0.0;;
702 aij = ai;
703 jx = kx;
704 for (j = 0; j < lenx; j++) {
705 a_elem[0] = a_i[aij];
706 a_elem[1] = a_i[aij + 1];
707
708 x_elem[0] = x_i[jx];
709 x_elem[1] = x_i[jx + 1];
710 {
711 prod[0] =
712 (double) a_elem[0] * x_elem[0] -
713 (double) a_elem[1] * x_elem[1];
714 prod[1] =
715 (double) a_elem[0] * x_elem[1] +
716 (double) a_elem[1] * x_elem[0];
717 }
718 sum[0] = sum[0] + prod[0];
719 sum[1] = sum[1] + prod[1];
720 aij += incaij;
721 jx += incx;
722 }
723 {
724 tmp1[0] =
725 (double) sum[0] * alpha_i[0] - (double) sum[1] * alpha_i[1];
726 tmp1[1] =
727 (double) sum[0] * alpha_i[1] + (double) sum[1] * alpha_i[0];
728 }
729 y_elem[0] = y_i[iy];
730 y_elem[1] = y_i[iy + 1];
731 {
732 tmp2[0] =
733 (double) y_elem[0] * beta_i[0] -
734 (double) y_elem[1] * beta_i[1];
735 tmp2[1] =
736 (double) y_elem[0] * beta_i[1] +
737 (double) y_elem[1] * beta_i[0];
738 }
739 tmp1[0] = tmp1[0] + tmp2[0];
740 tmp1[1] = tmp1[1] + tmp2[1];
741 y_i[iy] = tmp1[0];
742 y_i[iy + 1] = tmp1[1];
743 ai += incai;
744 iy += incy;
745 }
746 }
747
748 }
749 }
750
751
752
753 break;
754 }
755 case blas_prec_extra:{
756
757 int i, j;
758 int iy, jx, kx, ky;
759 int lenx, leny;
760 int ai, aij;
761 int incai, incaij;
762
763 const float *a_i = (float *) a;
764 const float *x_i = (float *) x;
765 float *y_i = (float *) y;
766 float *alpha_i = (float *) alpha;
767 float *beta_i = (float *) beta;
768 float a_elem[2];
769 float x_elem[2];
770 float y_elem[2];
771 double head_prod[2], tail_prod[2];
772 double head_sum[2], tail_sum[2];
773 double head_tmp1[2], tail_tmp1[2];
774 double head_tmp2[2], tail_tmp2[2];
775 FPU_FIX_DECL;
776
777 /* all error calls */
778 if (m < 0)
779 BLAS_error(routine_name, -3, m, 0);
780 else if (n <= 0)
781 BLAS_error(routine_name, -4, n, 0);
782 else if (incx == 0)
783 BLAS_error(routine_name, -9, incx, 0);
784 else if (incy == 0)
785 BLAS_error(routine_name, -12, incy, 0);
786
787 if ((order == blas_rowmajor) && (trans == blas_no_trans)) {
788 lenx = n;
789 leny = m;
790 incai = lda;
791 incaij = 1;
792 } else if ((order == blas_rowmajor) && (trans != blas_no_trans)) {
793 lenx = m;
794 leny = n;
795 incai = 1;
796 incaij = lda;
797 } else if ((order == blas_colmajor) && (trans == blas_no_trans)) {
798 lenx = n;
799 leny = m;
800 incai = 1;
801 incaij = lda;
802 } else { /* colmajor and blas_trans */
803 lenx = m;
804 leny = n;
805 incai = lda;
806 incaij = 1;
807 }
808 if ((order == blas_colmajor && lda < m) ||
809 (order == blas_rowmajor && lda < n))
810 BLAS_error(routine_name, -7, lda, NULL);
811
812 FPU_FIX_START;
813
814 incx *= 2;
815 incy *= 2;
816 incai *= 2;
817 incaij *= 2;
818
819 if (incx > 0)
820 kx = 0;
821 else
822 kx = (1 - lenx) * incx;
823 if (incy > 0)
824 ky = 0;
825 else
826 ky = (1 - leny) * incy;
827
828 /* No extra-precision needed for alpha = 0 */
829 if (alpha_i[0] == 0.0 && alpha_i[1] == 0.0) {
830 if (beta_i[0] == 0.0 && beta_i[1] == 0.0) {
831 iy = ky;
832 for (i = 0; i < leny; i++) {
833 y_i[iy] = 0.0;
834 y_i[iy + 1] = 0.0;
835 iy += incy;
836 }
837 } else if (!(beta_i[0] == 0.0 && beta_i[1] == 0.0)) {
838 iy = ky;
839 for (i = 0; i < leny; i++) {
840 y_elem[0] = y_i[iy];
841 y_elem[1] = y_i[iy + 1];
842 {
843 double head_e1, tail_e1;
844 double d1;
845 double d2;
846 /* Real part */
847 d1 = (double) y_elem[0] * beta_i[0];
848 d2 = (double) -y_elem[1] * beta_i[1];
849 {
850 /* Compute double-double = double + double. */
851 double e, t1, t2;
852
853 /* Knuth trick. */
854 t1 = d1 + d2;
855 e = t1 - d1;
856 t2 = ((d2 - e) + (d1 - (t1 - e)));
857
858 /* The result is t1 + t2, after normalization. */
859 head_e1 = t1 + t2;
860 tail_e1 = t2 - (head_e1 - t1);
861 }
862 head_tmp1[0] = head_e1;
863 tail_tmp1[0] = tail_e1;
864 /* imaginary part */
865 d1 = (double) y_elem[0] * beta_i[1];
866 d2 = (double) y_elem[1] * beta_i[0];
867 {
868 /* Compute double-double = double + double. */
869 double e, t1, t2;
870
871 /* Knuth trick. */
872 t1 = d1 + d2;
873 e = t1 - d1;
874 t2 = ((d2 - e) + (d1 - (t1 - e)));
875
876 /* The result is t1 + t2, after normalization. */
877 head_e1 = t1 + t2;
878 tail_e1 = t2 - (head_e1 - t1);
879 }
880 head_tmp1[1] = head_e1;
881 tail_tmp1[1] = tail_e1;
882 }
883 y_i[iy] = head_tmp1[0];
884 y_i[iy + 1] = head_tmp1[1];
885 iy += incy;
886 }
887 }
888 } else {
889 if (trans == blas_conj_trans) {
890
891 /* if beta = 0, we can save m multiplies: y = alpha*A*x */
892 if (beta_i[0] == 0.0 && beta_i[1] == 0.0) {
893 /* save m more multiplies if alpha = 1 */
894 if ((alpha_i[0] == 1.0 && alpha_i[1] == 0.0)) {
895 ai = 0;
896 iy = ky;
897 for (i = 0; i < leny; i++) {
898 head_sum[0] = head_sum[1] = tail_sum[0] = tail_sum[1] = 0.0;
899 aij = ai;
900 jx = kx;
901 for (j = 0; j < lenx; j++) {
902 a_elem[0] = a_i[aij];
903 a_elem[1] = a_i[aij + 1];
904 a_elem[1] = -a_elem[1];
905 x_elem[0] = x_i[jx];
906 x_elem[1] = x_i[jx + 1];
907 {
908 double head_e1, tail_e1;
909 double d1;
910 double d2;
911 /* Real part */
912 d1 = (double) a_elem[0] * x_elem[0];
913 d2 = (double) -a_elem[1] * x_elem[1];
914 {
915 /* Compute double-double = double + double. */
916 double e, t1, t2;
917
918 /* Knuth trick. */
919 t1 = d1 + d2;
920 e = t1 - d1;
921 t2 = ((d2 - e) + (d1 - (t1 - e)));
922
923 /* The result is t1 + t2, after normalization. */
924 head_e1 = t1 + t2;
925 tail_e1 = t2 - (head_e1 - t1);
926 }
927 head_prod[0] = head_e1;
928 tail_prod[0] = tail_e1;
929 /* imaginary part */
930 d1 = (double) a_elem[0] * x_elem[1];
931 d2 = (double) a_elem[1] * x_elem[0];
932 {
933 /* Compute double-double = double + double. */
934 double e, t1, t2;
935
936 /* Knuth trick. */
937 t1 = d1 + d2;
938 e = t1 - d1;
939 t2 = ((d2 - e) + (d1 - (t1 - e)));
940
941 /* The result is t1 + t2, after normalization. */
942 head_e1 = t1 + t2;
943 tail_e1 = t2 - (head_e1 - t1);
944 }
945 head_prod[1] = head_e1;
946 tail_prod[1] = tail_e1;
947 }
948 {
949 double head_t, tail_t;
950 double head_a, tail_a;
951 double head_b, tail_b;
952 /* Real part */
953 head_a = head_sum[0];
954 tail_a = tail_sum[0];
955 head_b = head_prod[0];
956 tail_b = tail_prod[0];
957 {
958 /* Compute double-double = double-double + double-double. */
959 double bv;
960 double s1, s2, t1, t2;
961
962 /* Add two hi words. */
963 s1 = head_a + head_b;
964 bv = s1 - head_a;
965 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
966
967 /* Add two lo words. */
968 t1 = tail_a + tail_b;
969 bv = t1 - tail_a;
970 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
971
972 s2 += t1;
973
974 /* Renormalize (s1, s2) to (t1, s2) */
975 t1 = s1 + s2;
976 s2 = s2 - (t1 - s1);
977
978 t2 += s2;
979
980 /* Renormalize (t1, t2) */
981 head_t = t1 + t2;
982 tail_t = t2 - (head_t - t1);
983 }
984 head_sum[0] = head_t;
985 tail_sum[0] = tail_t;
986 /* Imaginary part */
987 head_a = head_sum[1];
988 tail_a = tail_sum[1];
989 head_b = head_prod[1];
990 tail_b = tail_prod[1];
991 {
992 /* Compute double-double = double-double + double-double. */
993 double bv;
994 double s1, s2, t1, t2;
995
996 /* Add two hi words. */
997 s1 = head_a + head_b;
998 bv = s1 - head_a;
999 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1000
1001 /* Add two lo words. */
1002 t1 = tail_a + tail_b;
1003 bv = t1 - tail_a;
1004 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1005
1006 s2 += t1;
1007
1008 /* Renormalize (s1, s2) to (t1, s2) */
1009 t1 = s1 + s2;
1010 s2 = s2 - (t1 - s1);
1011
1012 t2 += s2;
1013
1014 /* Renormalize (t1, t2) */
1015 head_t = t1 + t2;
1016 tail_t = t2 - (head_t - t1);
1017 }
1018 head_sum[1] = head_t;
1019 tail_sum[1] = tail_t;
1020 }
1021 aij += incaij;
1022 jx += incx;
1023 }
1024 y_i[iy] = head_sum[0];
1025 y_i[iy + 1] = head_sum[1];
1026 ai += incai;
1027 iy += incy;
1028 }
1029 } else {
1030 ai = 0;
1031 iy = ky;
1032 for (i = 0; i < leny; i++) {
1033 head_sum[0] = head_sum[1] = tail_sum[0] = tail_sum[1] = 0.0;
1034 aij = ai;
1035 jx = kx;
1036 for (j = 0; j < lenx; j++) {
1037 a_elem[0] = a_i[aij];
1038 a_elem[1] = a_i[aij + 1];
1039 a_elem[1] = -a_elem[1];
1040 x_elem[0] = x_i[jx];
1041 x_elem[1] = x_i[jx + 1];
1042 {
1043 double head_e1, tail_e1;
1044 double d1;
1045 double d2;
1046 /* Real part */
1047 d1 = (double) a_elem[0] * x_elem[0];
1048 d2 = (double) -a_elem[1] * x_elem[1];
1049 {
1050 /* Compute double-double = double + double. */
1051 double e, t1, t2;
1052
1053 /* Knuth trick. */
1054 t1 = d1 + d2;
1055 e = t1 - d1;
1056 t2 = ((d2 - e) + (d1 - (t1 - e)));
1057
1058 /* The result is t1 + t2, after normalization. */
1059 head_e1 = t1 + t2;
1060 tail_e1 = t2 - (head_e1 - t1);
1061 }
1062 head_prod[0] = head_e1;
1063 tail_prod[0] = tail_e1;
1064 /* imaginary part */
1065 d1 = (double) a_elem[0] * x_elem[1];
1066 d2 = (double) a_elem[1] * x_elem[0];
1067 {
1068 /* Compute double-double = double + double. */
1069 double e, t1, t2;
1070
1071 /* Knuth trick. */
1072 t1 = d1 + d2;
1073 e = t1 - d1;
1074 t2 = ((d2 - e) + (d1 - (t1 - e)));
1075
1076 /* The result is t1 + t2, after normalization. */
1077 head_e1 = t1 + t2;
1078 tail_e1 = t2 - (head_e1 - t1);
1079 }
1080 head_prod[1] = head_e1;
1081 tail_prod[1] = tail_e1;
1082 }
1083 {
1084 double head_t, tail_t;
1085 double head_a, tail_a;
1086 double head_b, tail_b;
1087 /* Real part */
1088 head_a = head_sum[0];
1089 tail_a = tail_sum[0];
1090 head_b = head_prod[0];
1091 tail_b = tail_prod[0];
1092 {
1093 /* Compute double-double = double-double + double-double. */
1094 double bv;
1095 double s1, s2, t1, t2;
1096
1097 /* Add two hi words. */
1098 s1 = head_a + head_b;
1099 bv = s1 - head_a;
1100 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1101
1102 /* Add two lo words. */
1103 t1 = tail_a + tail_b;
1104 bv = t1 - tail_a;
1105 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1106
1107 s2 += t1;
1108
1109 /* Renormalize (s1, s2) to (t1, s2) */
1110 t1 = s1 + s2;
1111 s2 = s2 - (t1 - s1);
1112
1113 t2 += s2;
1114
1115 /* Renormalize (t1, t2) */
1116 head_t = t1 + t2;
1117 tail_t = t2 - (head_t - t1);
1118 }
1119 head_sum[0] = head_t;
1120 tail_sum[0] = tail_t;
1121 /* Imaginary part */
1122 head_a = head_sum[1];
1123 tail_a = tail_sum[1];
1124 head_b = head_prod[1];
1125 tail_b = tail_prod[1];
1126 {
1127 /* Compute double-double = double-double + double-double. */
1128 double bv;
1129 double s1, s2, t1, t2;
1130
1131 /* Add two hi words. */
1132 s1 = head_a + head_b;
1133 bv = s1 - head_a;
1134 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1135
1136 /* Add two lo words. */
1137 t1 = tail_a + tail_b;
1138 bv = t1 - tail_a;
1139 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1140
1141 s2 += t1;
1142
1143 /* Renormalize (s1, s2) to (t1, s2) */
1144 t1 = s1 + s2;
1145 s2 = s2 - (t1 - s1);
1146
1147 t2 += s2;
1148
1149 /* Renormalize (t1, t2) */
1150 head_t = t1 + t2;
1151 tail_t = t2 - (head_t - t1);
1152 }
1153 head_sum[1] = head_t;
1154 tail_sum[1] = tail_t;
1155 }
1156 aij += incaij;
1157 jx += incx;
1158 }
1159 {
1160 double cd[2];
1161 cd[0] = (double) alpha_i[0];
1162 cd[1] = (double) alpha_i[1];
1163 {
1164 /* Compute complex-extra = complex-extra * complex-double. */
1165 double head_a0, tail_a0;
1166 double head_a1, tail_a1;
1167 double head_t1, tail_t1;
1168 double head_t2, tail_t2;
1169 head_a0 = head_sum[0];
1170 tail_a0 = tail_sum[0];
1171 head_a1 = head_sum[1];
1172 tail_a1 = tail_sum[1];
1173 /* real part */
1174 {
1175 /* Compute double-double = double-double * double. */
1176 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
1177
1178 con = head_a0 * split;
1179 a11 = con - head_a0;
1180 a11 = con - a11;
1181 a21 = head_a0 - a11;
1182 con = cd[0] * split;
1183 b1 = con - cd[0];
1184 b1 = con - b1;
1185 b2 = cd[0] - b1;
1186
1187 c11 = head_a0 * cd[0];
1188 c21 =
1189 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
1190
1191 c2 = tail_a0 * cd[0];
1192 t1 = c11 + c2;
1193 t2 = (c2 - (t1 - c11)) + c21;
1194
1195 head_t1 = t1 + t2;
1196 tail_t1 = t2 - (head_t1 - t1);
1197 }
1198 {
1199 /* Compute double-double = double-double * double. */
1200 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
1201
1202 con = head_a1 * split;
1203 a11 = con - head_a1;
1204 a11 = con - a11;
1205 a21 = head_a1 - a11;
1206 con = cd[1] * split;
1207 b1 = con - cd[1];
1208 b1 = con - b1;
1209 b2 = cd[1] - b1;
1210
1211 c11 = head_a1 * cd[1];
1212 c21 =
1213 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
1214
1215 c2 = tail_a1 * cd[1];
1216 t1 = c11 + c2;
1217 t2 = (c2 - (t1 - c11)) + c21;
1218
1219 head_t2 = t1 + t2;
1220 tail_t2 = t2 - (head_t2 - t1);
1221 }
1222 head_t2 = -head_t2;
1223 tail_t2 = -tail_t2;
1224 {
1225 /* Compute double-double = double-double + double-double. */
1226 double bv;
1227 double s1, s2, t1, t2;
1228
1229 /* Add two hi words. */
1230 s1 = head_t1 + head_t2;
1231 bv = s1 - head_t1;
1232 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
1233
1234 /* Add two lo words. */
1235 t1 = tail_t1 + tail_t2;
1236 bv = t1 - tail_t1;
1237 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
1238
1239 s2 += t1;
1240
1241 /* Renormalize (s1, s2) to (t1, s2) */
1242 t1 = s1 + s2;
1243 s2 = s2 - (t1 - s1);
1244
1245 t2 += s2;
1246
1247 /* Renormalize (t1, t2) */
1248 head_t1 = t1 + t2;
1249 tail_t1 = t2 - (head_t1 - t1);
1250 }
1251 head_tmp1[0] = head_t1;
1252 tail_tmp1[0] = tail_t1;
1253 /* imaginary part */
1254 {
1255 /* Compute double-double = double-double * double. */
1256 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
1257
1258 con = head_a1 * split;
1259 a11 = con - head_a1;
1260 a11 = con - a11;
1261 a21 = head_a1 - a11;
1262 con = cd[0] * split;
1263 b1 = con - cd[0];
1264 b1 = con - b1;
1265 b2 = cd[0] - b1;
1266
1267 c11 = head_a1 * cd[0];
1268 c21 =
1269 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
1270
1271 c2 = tail_a1 * cd[0];
1272 t1 = c11 + c2;
1273 t2 = (c2 - (t1 - c11)) + c21;
1274
1275 head_t1 = t1 + t2;
1276 tail_t1 = t2 - (head_t1 - t1);
1277 }
1278 {
1279 /* Compute double-double = double-double * double. */
1280 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
1281
1282 con = head_a0 * split;
1283 a11 = con - head_a0;
1284 a11 = con - a11;
1285 a21 = head_a0 - a11;
1286 con = cd[1] * split;
1287 b1 = con - cd[1];
1288 b1 = con - b1;
1289 b2 = cd[1] - b1;
1290
1291 c11 = head_a0 * cd[1];
1292 c21 =
1293 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
1294
1295 c2 = tail_a0 * cd[1];
1296 t1 = c11 + c2;
1297 t2 = (c2 - (t1 - c11)) + c21;
1298
1299 head_t2 = t1 + t2;
1300 tail_t2 = t2 - (head_t2 - t1);
1301 }
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_t1 + head_t2;
1309 bv = s1 - head_t1;
1310 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
1311
1312 /* Add two lo words. */
1313 t1 = tail_t1 + tail_t2;
1314 bv = t1 - tail_t1;
1315 t2 = ((tail_t2 - bv) + (tail_t1 - (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_t1 = t1 + t2;
1327 tail_t1 = t2 - (head_t1 - t1);
1328 }
1329 head_tmp1[1] = head_t1;
1330 tail_tmp1[1] = tail_t1;
1331 }
1332
1333 }
1334 y_i[iy] = head_tmp1[0];
1335 y_i[iy + 1] = head_tmp1[1];
1336 ai += incai;
1337 iy += incy;
1338 }
1339 }
1340 } else {
1341 /* the most general form, y = alpha*A*x + beta*y */
1342 ai = 0;
1343 iy = ky;
1344 for (i = 0; i < leny; i++) {
1345 head_sum[0] = head_sum[1] = tail_sum[0] = tail_sum[1] = 0.0;;
1346 aij = ai;
1347 jx = kx;
1348 for (j = 0; j < lenx; j++) {
1349 a_elem[0] = a_i[aij];
1350 a_elem[1] = a_i[aij + 1];
1351 a_elem[1] = -a_elem[1];
1352 x_elem[0] = x_i[jx];
1353 x_elem[1] = x_i[jx + 1];
1354 {
1355 double head_e1, tail_e1;
1356 double d1;
1357 double d2;
1358 /* Real part */
1359 d1 = (double) a_elem[0] * x_elem[0];
1360 d2 = (double) -a_elem[1] * x_elem[1];
1361 {
1362 /* Compute double-double = double + double. */
1363 double e, t1, t2;
1364
1365 /* Knuth trick. */
1366 t1 = d1 + d2;
1367 e = t1 - d1;
1368 t2 = ((d2 - e) + (d1 - (t1 - e)));
1369
1370 /* The result is t1 + t2, after normalization. */
1371 head_e1 = t1 + t2;
1372 tail_e1 = t2 - (head_e1 - t1);
1373 }
1374 head_prod[0] = head_e1;
1375 tail_prod[0] = tail_e1;
1376 /* imaginary part */
1377 d1 = (double) a_elem[0] * x_elem[1];
1378 d2 = (double) a_elem[1] * x_elem[0];
1379 {
1380 /* Compute double-double = double + double. */
1381 double e, t1, t2;
1382
1383 /* Knuth trick. */
1384 t1 = d1 + d2;
1385 e = t1 - d1;
1386 t2 = ((d2 - e) + (d1 - (t1 - e)));
1387
1388 /* The result is t1 + t2, after normalization. */
1389 head_e1 = t1 + t2;
1390 tail_e1 = t2 - (head_e1 - t1);
1391 }
1392 head_prod[1] = head_e1;
1393 tail_prod[1] = tail_e1;
1394 }
1395 {
1396 double head_t, tail_t;
1397 double head_a, tail_a;
1398 double head_b, tail_b;
1399 /* Real part */
1400 head_a = head_sum[0];
1401 tail_a = tail_sum[0];
1402 head_b = head_prod[0];
1403 tail_b = tail_prod[0];
1404 {
1405 /* Compute double-double = double-double + double-double. */
1406 double bv;
1407 double s1, s2, t1, t2;
1408
1409 /* Add two hi words. */
1410 s1 = head_a + head_b;
1411 bv = s1 - head_a;
1412 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1413
1414 /* Add two lo words. */
1415 t1 = tail_a + tail_b;
1416 bv = t1 - tail_a;
1417 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1418
1419 s2 += t1;
1420
1421 /* Renormalize (s1, s2) to (t1, s2) */
1422 t1 = s1 + s2;
1423 s2 = s2 - (t1 - s1);
1424
1425 t2 += s2;
1426
1427 /* Renormalize (t1, t2) */
1428 head_t = t1 + t2;
1429 tail_t = t2 - (head_t - t1);
1430 }
1431 head_sum[0] = head_t;
1432 tail_sum[0] = tail_t;
1433 /* Imaginary part */
1434 head_a = head_sum[1];
1435 tail_a = tail_sum[1];
1436 head_b = head_prod[1];
1437 tail_b = tail_prod[1];
1438 {
1439 /* Compute double-double = double-double + double-double. */
1440 double bv;
1441 double s1, s2, t1, t2;
1442
1443 /* Add two hi words. */
1444 s1 = head_a + head_b;
1445 bv = s1 - head_a;
1446 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1447
1448 /* Add two lo words. */
1449 t1 = tail_a + tail_b;
1450 bv = t1 - tail_a;
1451 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1452
1453 s2 += t1;
1454
1455 /* Renormalize (s1, s2) to (t1, s2) */
1456 t1 = s1 + s2;
1457 s2 = s2 - (t1 - s1);
1458
1459 t2 += s2;
1460
1461 /* Renormalize (t1, t2) */
1462 head_t = t1 + t2;
1463 tail_t = t2 - (head_t - t1);
1464 }
1465 head_sum[1] = head_t;
1466 tail_sum[1] = tail_t;
1467 }
1468 aij += incaij;
1469 jx += incx;
1470 }
1471 {
1472 double cd[2];
1473 cd[0] = (double) alpha_i[0];
1474 cd[1] = (double) alpha_i[1];
1475 {
1476 /* Compute complex-extra = complex-extra * complex-double. */
1477 double head_a0, tail_a0;
1478 double head_a1, tail_a1;
1479 double head_t1, tail_t1;
1480 double head_t2, tail_t2;
1481 head_a0 = head_sum[0];
1482 tail_a0 = tail_sum[0];
1483 head_a1 = head_sum[1];
1484 tail_a1 = tail_sum[1];
1485 /* real part */
1486 {
1487 /* Compute double-double = double-double * double. */
1488 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
1489
1490 con = head_a0 * split;
1491 a11 = con - head_a0;
1492 a11 = con - a11;
1493 a21 = head_a0 - a11;
1494 con = cd[0] * split;
1495 b1 = con - cd[0];
1496 b1 = con - b1;
1497 b2 = cd[0] - b1;
1498
1499 c11 = head_a0 * cd[0];
1500 c21 =
1501 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
1502
1503 c2 = tail_a0 * cd[0];
1504 t1 = c11 + c2;
1505 t2 = (c2 - (t1 - c11)) + c21;
1506
1507 head_t1 = t1 + t2;
1508 tail_t1 = t2 - (head_t1 - t1);
1509 }
1510 {
1511 /* Compute double-double = double-double * double. */
1512 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
1513
1514 con = head_a1 * split;
1515 a11 = con - head_a1;
1516 a11 = con - a11;
1517 a21 = head_a1 - a11;
1518 con = cd[1] * split;
1519 b1 = con - cd[1];
1520 b1 = con - b1;
1521 b2 = cd[1] - b1;
1522
1523 c11 = head_a1 * cd[1];
1524 c21 =
1525 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
1526
1527 c2 = tail_a1 * cd[1];
1528 t1 = c11 + c2;
1529 t2 = (c2 - (t1 - c11)) + c21;
1530
1531 head_t2 = t1 + t2;
1532 tail_t2 = t2 - (head_t2 - t1);
1533 }
1534 head_t2 = -head_t2;
1535 tail_t2 = -tail_t2;
1536 {
1537 /* Compute double-double = double-double + double-double. */
1538 double bv;
1539 double s1, s2, t1, t2;
1540
1541 /* Add two hi words. */
1542 s1 = head_t1 + head_t2;
1543 bv = s1 - head_t1;
1544 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
1545
1546 /* Add two lo words. */
1547 t1 = tail_t1 + tail_t2;
1548 bv = t1 - tail_t1;
1549 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
1550
1551 s2 += t1;
1552
1553 /* Renormalize (s1, s2) to (t1, s2) */
1554 t1 = s1 + s2;
1555 s2 = s2 - (t1 - s1);
1556
1557 t2 += s2;
1558
1559 /* Renormalize (t1, t2) */
1560 head_t1 = t1 + t2;
1561 tail_t1 = t2 - (head_t1 - t1);
1562 }
1563 head_tmp1[0] = head_t1;
1564 tail_tmp1[0] = tail_t1;
1565 /* imaginary part */
1566 {
1567 /* Compute double-double = double-double * double. */
1568 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
1569
1570 con = head_a1 * split;
1571 a11 = con - head_a1;
1572 a11 = con - a11;
1573 a21 = head_a1 - a11;
1574 con = cd[0] * split;
1575 b1 = con - cd[0];
1576 b1 = con - b1;
1577 b2 = cd[0] - b1;
1578
1579 c11 = head_a1 * cd[0];
1580 c21 =
1581 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
1582
1583 c2 = tail_a1 * cd[0];
1584 t1 = c11 + c2;
1585 t2 = (c2 - (t1 - c11)) + c21;
1586
1587 head_t1 = t1 + t2;
1588 tail_t1 = t2 - (head_t1 - t1);
1589 }
1590 {
1591 /* Compute double-double = double-double * double. */
1592 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
1593
1594 con = head_a0 * split;
1595 a11 = con - head_a0;
1596 a11 = con - a11;
1597 a21 = head_a0 - a11;
1598 con = cd[1] * split;
1599 b1 = con - cd[1];
1600 b1 = con - b1;
1601 b2 = cd[1] - b1;
1602
1603 c11 = head_a0 * cd[1];
1604 c21 =
1605 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
1606
1607 c2 = tail_a0 * cd[1];
1608 t1 = c11 + c2;
1609 t2 = (c2 - (t1 - c11)) + c21;
1610
1611 head_t2 = t1 + t2;
1612 tail_t2 = t2 - (head_t2 - t1);
1613 }
1614 {
1615 /* Compute double-double = double-double + double-double. */
1616 double bv;
1617 double s1, s2, t1, t2;
1618
1619 /* Add two hi words. */
1620 s1 = head_t1 + head_t2;
1621 bv = s1 - head_t1;
1622 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
1623
1624 /* Add two lo words. */
1625 t1 = tail_t1 + tail_t2;
1626 bv = t1 - tail_t1;
1627 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
1628
1629 s2 += t1;
1630
1631 /* Renormalize (s1, s2) to (t1, s2) */
1632 t1 = s1 + s2;
1633 s2 = s2 - (t1 - s1);
1634
1635 t2 += s2;
1636
1637 /* Renormalize (t1, t2) */
1638 head_t1 = t1 + t2;
1639 tail_t1 = t2 - (head_t1 - t1);
1640 }
1641 head_tmp1[1] = head_t1;
1642 tail_tmp1[1] = tail_t1;
1643 }
1644
1645 }
1646 y_elem[0] = y_i[iy];
1647 y_elem[1] = y_i[iy + 1];
1648 {
1649 double head_e1, tail_e1;
1650 double d1;
1651 double d2;
1652 /* Real part */
1653 d1 = (double) y_elem[0] * beta_i[0];
1654 d2 = (double) -y_elem[1] * beta_i[1];
1655 {
1656 /* Compute double-double = double + double. */
1657 double e, t1, t2;
1658
1659 /* Knuth trick. */
1660 t1 = d1 + d2;
1661 e = t1 - d1;
1662 t2 = ((d2 - e) + (d1 - (t1 - e)));
1663
1664 /* The result is t1 + t2, after normalization. */
1665 head_e1 = t1 + t2;
1666 tail_e1 = t2 - (head_e1 - t1);
1667 }
1668 head_tmp2[0] = head_e1;
1669 tail_tmp2[0] = tail_e1;
1670 /* imaginary part */
1671 d1 = (double) y_elem[0] * beta_i[1];
1672 d2 = (double) y_elem[1] * beta_i[0];
1673 {
1674 /* Compute double-double = double + double. */
1675 double e, t1, t2;
1676
1677 /* Knuth trick. */
1678 t1 = d1 + d2;
1679 e = t1 - d1;
1680 t2 = ((d2 - e) + (d1 - (t1 - e)));
1681
1682 /* The result is t1 + t2, after normalization. */
1683 head_e1 = t1 + t2;
1684 tail_e1 = t2 - (head_e1 - t1);
1685 }
1686 head_tmp2[1] = head_e1;
1687 tail_tmp2[1] = tail_e1;
1688 }
1689 {
1690 double head_t, tail_t;
1691 double head_a, tail_a;
1692 double head_b, tail_b;
1693 /* Real part */
1694 head_a = head_tmp1[0];
1695 tail_a = tail_tmp1[0];
1696 head_b = head_tmp2[0];
1697 tail_b = tail_tmp2[0];
1698 {
1699 /* Compute double-double = double-double + double-double. */
1700 double bv;
1701 double s1, s2, t1, t2;
1702
1703 /* Add two hi words. */
1704 s1 = head_a + head_b;
1705 bv = s1 - head_a;
1706 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1707
1708 /* Add two lo words. */
1709 t1 = tail_a + tail_b;
1710 bv = t1 - tail_a;
1711 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1712
1713 s2 += t1;
1714
1715 /* Renormalize (s1, s2) to (t1, s2) */
1716 t1 = s1 + s2;
1717 s2 = s2 - (t1 - s1);
1718
1719 t2 += s2;
1720
1721 /* Renormalize (t1, t2) */
1722 head_t = t1 + t2;
1723 tail_t = t2 - (head_t - t1);
1724 }
1725 head_tmp1[0] = head_t;
1726 tail_tmp1[0] = tail_t;
1727 /* Imaginary part */
1728 head_a = head_tmp1[1];
1729 tail_a = tail_tmp1[1];
1730 head_b = head_tmp2[1];
1731 tail_b = tail_tmp2[1];
1732 {
1733 /* Compute double-double = double-double + double-double. */
1734 double bv;
1735 double s1, s2, t1, t2;
1736
1737 /* Add two hi words. */
1738 s1 = head_a + head_b;
1739 bv = s1 - head_a;
1740 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1741
1742 /* Add two lo words. */
1743 t1 = tail_a + tail_b;
1744 bv = t1 - tail_a;
1745 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1746
1747 s2 += t1;
1748
1749 /* Renormalize (s1, s2) to (t1, s2) */
1750 t1 = s1 + s2;
1751 s2 = s2 - (t1 - s1);
1752
1753 t2 += s2;
1754
1755 /* Renormalize (t1, t2) */
1756 head_t = t1 + t2;
1757 tail_t = t2 - (head_t - t1);
1758 }
1759 head_tmp1[1] = head_t;
1760 tail_tmp1[1] = tail_t;
1761 }
1762 y_i[iy] = head_tmp1[0];
1763 y_i[iy + 1] = head_tmp1[1];
1764 ai += incai;
1765 iy += incy;
1766 }
1767 }
1768
1769 } else {
1770
1771 /* if beta = 0, we can save m multiplies: y = alpha*A*x */
1772 if (beta_i[0] == 0.0 && beta_i[1] == 0.0) {
1773 /* save m more multiplies if alpha = 1 */
1774 if ((alpha_i[0] == 1.0 && alpha_i[1] == 0.0)) {
1775 ai = 0;
1776 iy = ky;
1777 for (i = 0; i < leny; i++) {
1778 head_sum[0] = head_sum[1] = tail_sum[0] = tail_sum[1] = 0.0;
1779 aij = ai;
1780 jx = kx;
1781 for (j = 0; j < lenx; j++) {
1782 a_elem[0] = a_i[aij];
1783 a_elem[1] = a_i[aij + 1];
1784
1785 x_elem[0] = x_i[jx];
1786 x_elem[1] = x_i[jx + 1];
1787 {
1788 double head_e1, tail_e1;
1789 double d1;
1790 double d2;
1791 /* Real part */
1792 d1 = (double) a_elem[0] * x_elem[0];
1793 d2 = (double) -a_elem[1] * x_elem[1];
1794 {
1795 /* Compute double-double = double + double. */
1796 double e, t1, t2;
1797
1798 /* Knuth trick. */
1799 t1 = d1 + d2;
1800 e = t1 - d1;
1801 t2 = ((d2 - e) + (d1 - (t1 - e)));
1802
1803 /* The result is t1 + t2, after normalization. */
1804 head_e1 = t1 + t2;
1805 tail_e1 = t2 - (head_e1 - t1);
1806 }
1807 head_prod[0] = head_e1;
1808 tail_prod[0] = tail_e1;
1809 /* imaginary part */
1810 d1 = (double) a_elem[0] * x_elem[1];
1811 d2 = (double) a_elem[1] * x_elem[0];
1812 {
1813 /* Compute double-double = double + double. */
1814 double e, t1, t2;
1815
1816 /* Knuth trick. */
1817 t1 = d1 + d2;
1818 e = t1 - d1;
1819 t2 = ((d2 - e) + (d1 - (t1 - e)));
1820
1821 /* The result is t1 + t2, after normalization. */
1822 head_e1 = t1 + t2;
1823 tail_e1 = t2 - (head_e1 - t1);
1824 }
1825 head_prod[1] = head_e1;
1826 tail_prod[1] = tail_e1;
1827 }
1828 {
1829 double head_t, tail_t;
1830 double head_a, tail_a;
1831 double head_b, tail_b;
1832 /* Real part */
1833 head_a = head_sum[0];
1834 tail_a = tail_sum[0];
1835 head_b = head_prod[0];
1836 tail_b = tail_prod[0];
1837 {
1838 /* Compute double-double = double-double + double-double. */
1839 double bv;
1840 double s1, s2, t1, t2;
1841
1842 /* Add two hi words. */
1843 s1 = head_a + head_b;
1844 bv = s1 - head_a;
1845 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1846
1847 /* Add two lo words. */
1848 t1 = tail_a + tail_b;
1849 bv = t1 - tail_a;
1850 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1851
1852 s2 += t1;
1853
1854 /* Renormalize (s1, s2) to (t1, s2) */
1855 t1 = s1 + s2;
1856 s2 = s2 - (t1 - s1);
1857
1858 t2 += s2;
1859
1860 /* Renormalize (t1, t2) */
1861 head_t = t1 + t2;
1862 tail_t = t2 - (head_t - t1);
1863 }
1864 head_sum[0] = head_t;
1865 tail_sum[0] = tail_t;
1866 /* Imaginary part */
1867 head_a = head_sum[1];
1868 tail_a = tail_sum[1];
1869 head_b = head_prod[1];
1870 tail_b = tail_prod[1];
1871 {
1872 /* Compute double-double = double-double + double-double. */
1873 double bv;
1874 double s1, s2, t1, t2;
1875
1876 /* Add two hi words. */
1877 s1 = head_a + head_b;
1878 bv = s1 - head_a;
1879 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1880
1881 /* Add two lo words. */
1882 t1 = tail_a + tail_b;
1883 bv = t1 - tail_a;
1884 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1885
1886 s2 += t1;
1887
1888 /* Renormalize (s1, s2) to (t1, s2) */
1889 t1 = s1 + s2;
1890 s2 = s2 - (t1 - s1);
1891
1892 t2 += s2;
1893
1894 /* Renormalize (t1, t2) */
1895 head_t = t1 + t2;
1896 tail_t = t2 - (head_t - t1);
1897 }
1898 head_sum[1] = head_t;
1899 tail_sum[1] = tail_t;
1900 }
1901 aij += incaij;
1902 jx += incx;
1903 }
1904 y_i[iy] = head_sum[0];
1905 y_i[iy + 1] = head_sum[1];
1906 ai += incai;
1907 iy += incy;
1908 }
1909 } else {
1910 ai = 0;
1911 iy = ky;
1912 for (i = 0; i < leny; i++) {
1913 head_sum[0] = head_sum[1] = tail_sum[0] = tail_sum[1] = 0.0;
1914 aij = ai;
1915 jx = kx;
1916 for (j = 0; j < lenx; j++) {
1917 a_elem[0] = a_i[aij];
1918 a_elem[1] = a_i[aij + 1];
1919
1920 x_elem[0] = x_i[jx];
1921 x_elem[1] = x_i[jx + 1];
1922 {
1923 double head_e1, tail_e1;
1924 double d1;
1925 double d2;
1926 /* Real part */
1927 d1 = (double) a_elem[0] * x_elem[0];
1928 d2 = (double) -a_elem[1] * x_elem[1];
1929 {
1930 /* Compute double-double = double + double. */
1931 double e, t1, t2;
1932
1933 /* Knuth trick. */
1934 t1 = d1 + d2;
1935 e = t1 - d1;
1936 t2 = ((d2 - e) + (d1 - (t1 - e)));
1937
1938 /* The result is t1 + t2, after normalization. */
1939 head_e1 = t1 + t2;
1940 tail_e1 = t2 - (head_e1 - t1);
1941 }
1942 head_prod[0] = head_e1;
1943 tail_prod[0] = tail_e1;
1944 /* imaginary part */
1945 d1 = (double) a_elem[0] * x_elem[1];
1946 d2 = (double) a_elem[1] * x_elem[0];
1947 {
1948 /* Compute double-double = double + double. */
1949 double e, t1, t2;
1950
1951 /* Knuth trick. */
1952 t1 = d1 + d2;
1953 e = t1 - d1;
1954 t2 = ((d2 - e) + (d1 - (t1 - e)));
1955
1956 /* The result is t1 + t2, after normalization. */
1957 head_e1 = t1 + t2;
1958 tail_e1 = t2 - (head_e1 - t1);
1959 }
1960 head_prod[1] = head_e1;
1961 tail_prod[1] = tail_e1;
1962 }
1963 {
1964 double head_t, tail_t;
1965 double head_a, tail_a;
1966 double head_b, tail_b;
1967 /* Real part */
1968 head_a = head_sum[0];
1969 tail_a = tail_sum[0];
1970 head_b = head_prod[0];
1971 tail_b = tail_prod[0];
1972 {
1973 /* Compute double-double = double-double + double-double. */
1974 double bv;
1975 double s1, s2, t1, t2;
1976
1977 /* Add two hi words. */
1978 s1 = head_a + head_b;
1979 bv = s1 - head_a;
1980 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1981
1982 /* Add two lo words. */
1983 t1 = tail_a + tail_b;
1984 bv = t1 - tail_a;
1985 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1986
1987 s2 += t1;
1988
1989 /* Renormalize (s1, s2) to (t1, s2) */
1990 t1 = s1 + s2;
1991 s2 = s2 - (t1 - s1);
1992
1993 t2 += s2;
1994
1995 /* Renormalize (t1, t2) */
1996 head_t = t1 + t2;
1997 tail_t = t2 - (head_t - t1);
1998 }
1999 head_sum[0] = head_t;
2000 tail_sum[0] = tail_t;
2001 /* Imaginary part */
2002 head_a = head_sum[1];
2003 tail_a = tail_sum[1];
2004 head_b = head_prod[1];
2005 tail_b = tail_prod[1];
2006 {
2007 /* Compute double-double = double-double + double-double. */
2008 double bv;
2009 double s1, s2, t1, t2;
2010
2011 /* Add two hi words. */
2012 s1 = head_a + head_b;
2013 bv = s1 - head_a;
2014 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2015
2016 /* Add two lo words. */
2017 t1 = tail_a + tail_b;
2018 bv = t1 - tail_a;
2019 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2020
2021 s2 += t1;
2022
2023 /* Renormalize (s1, s2) to (t1, s2) */
2024 t1 = s1 + s2;
2025 s2 = s2 - (t1 - s1);
2026
2027 t2 += s2;
2028
2029 /* Renormalize (t1, t2) */
2030 head_t = t1 + t2;
2031 tail_t = t2 - (head_t - t1);
2032 }
2033 head_sum[1] = head_t;
2034 tail_sum[1] = tail_t;
2035 }
2036 aij += incaij;
2037 jx += incx;
2038 }
2039 {
2040 double cd[2];
2041 cd[0] = (double) alpha_i[0];
2042 cd[1] = (double) alpha_i[1];
2043 {
2044 /* Compute complex-extra = complex-extra * complex-double. */
2045 double head_a0, tail_a0;
2046 double head_a1, tail_a1;
2047 double head_t1, tail_t1;
2048 double head_t2, tail_t2;
2049 head_a0 = head_sum[0];
2050 tail_a0 = tail_sum[0];
2051 head_a1 = head_sum[1];
2052 tail_a1 = tail_sum[1];
2053 /* real part */
2054 {
2055 /* Compute double-double = double-double * double. */
2056 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
2057
2058 con = head_a0 * split;
2059 a11 = con - head_a0;
2060 a11 = con - a11;
2061 a21 = head_a0 - a11;
2062 con = cd[0] * split;
2063 b1 = con - cd[0];
2064 b1 = con - b1;
2065 b2 = cd[0] - b1;
2066
2067 c11 = head_a0 * cd[0];
2068 c21 =
2069 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
2070
2071 c2 = tail_a0 * cd[0];
2072 t1 = c11 + c2;
2073 t2 = (c2 - (t1 - c11)) + c21;
2074
2075 head_t1 = t1 + t2;
2076 tail_t1 = t2 - (head_t1 - t1);
2077 }
2078 {
2079 /* Compute double-double = double-double * double. */
2080 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
2081
2082 con = head_a1 * split;
2083 a11 = con - head_a1;
2084 a11 = con - a11;
2085 a21 = head_a1 - a11;
2086 con = cd[1] * split;
2087 b1 = con - cd[1];
2088 b1 = con - b1;
2089 b2 = cd[1] - b1;
2090
2091 c11 = head_a1 * cd[1];
2092 c21 =
2093 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
2094
2095 c2 = tail_a1 * cd[1];
2096 t1 = c11 + c2;
2097 t2 = (c2 - (t1 - c11)) + c21;
2098
2099 head_t2 = t1 + t2;
2100 tail_t2 = t2 - (head_t2 - t1);
2101 }
2102 head_t2 = -head_t2;
2103 tail_t2 = -tail_t2;
2104 {
2105 /* Compute double-double = double-double + double-double. */
2106 double bv;
2107 double s1, s2, t1, t2;
2108
2109 /* Add two hi words. */
2110 s1 = head_t1 + head_t2;
2111 bv = s1 - head_t1;
2112 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
2113
2114 /* Add two lo words. */
2115 t1 = tail_t1 + tail_t2;
2116 bv = t1 - tail_t1;
2117 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
2118
2119 s2 += t1;
2120
2121 /* Renormalize (s1, s2) to (t1, s2) */
2122 t1 = s1 + s2;
2123 s2 = s2 - (t1 - s1);
2124
2125 t2 += s2;
2126
2127 /* Renormalize (t1, t2) */
2128 head_t1 = t1 + t2;
2129 tail_t1 = t2 - (head_t1 - t1);
2130 }
2131 head_tmp1[0] = head_t1;
2132 tail_tmp1[0] = tail_t1;
2133 /* imaginary part */
2134 {
2135 /* Compute double-double = double-double * double. */
2136 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
2137
2138 con = head_a1 * split;
2139 a11 = con - head_a1;
2140 a11 = con - a11;
2141 a21 = head_a1 - a11;
2142 con = cd[0] * split;
2143 b1 = con - cd[0];
2144 b1 = con - b1;
2145 b2 = cd[0] - b1;
2146
2147 c11 = head_a1 * cd[0];
2148 c21 =
2149 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
2150
2151 c2 = tail_a1 * cd[0];
2152 t1 = c11 + c2;
2153 t2 = (c2 - (t1 - c11)) + c21;
2154
2155 head_t1 = t1 + t2;
2156 tail_t1 = t2 - (head_t1 - t1);
2157 }
2158 {
2159 /* Compute double-double = double-double * double. */
2160 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
2161
2162 con = head_a0 * split;
2163 a11 = con - head_a0;
2164 a11 = con - a11;
2165 a21 = head_a0 - a11;
2166 con = cd[1] * split;
2167 b1 = con - cd[1];
2168 b1 = con - b1;
2169 b2 = cd[1] - b1;
2170
2171 c11 = head_a0 * cd[1];
2172 c21 =
2173 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
2174
2175 c2 = tail_a0 * cd[1];
2176 t1 = c11 + c2;
2177 t2 = (c2 - (t1 - c11)) + c21;
2178
2179 head_t2 = t1 + t2;
2180 tail_t2 = t2 - (head_t2 - t1);
2181 }
2182 {
2183 /* Compute double-double = double-double + double-double. */
2184 double bv;
2185 double s1, s2, t1, t2;
2186
2187 /* Add two hi words. */
2188 s1 = head_t1 + head_t2;
2189 bv = s1 - head_t1;
2190 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
2191
2192 /* Add two lo words. */
2193 t1 = tail_t1 + tail_t2;
2194 bv = t1 - tail_t1;
2195 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
2196
2197 s2 += t1;
2198
2199 /* Renormalize (s1, s2) to (t1, s2) */
2200 t1 = s1 + s2;
2201 s2 = s2 - (t1 - s1);
2202
2203 t2 += s2;
2204
2205 /* Renormalize (t1, t2) */
2206 head_t1 = t1 + t2;
2207 tail_t1 = t2 - (head_t1 - t1);
2208 }
2209 head_tmp1[1] = head_t1;
2210 tail_tmp1[1] = tail_t1;
2211 }
2212
2213 }
2214 y_i[iy] = head_tmp1[0];
2215 y_i[iy + 1] = head_tmp1[1];
2216 ai += incai;
2217 iy += incy;
2218 }
2219 }
2220 } else {
2221 /* the most general form, y = alpha*A*x + beta*y */
2222 ai = 0;
2223 iy = ky;
2224 for (i = 0; i < leny; i++) {
2225 head_sum[0] = head_sum[1] = tail_sum[0] = tail_sum[1] = 0.0;;
2226 aij = ai;
2227 jx = kx;
2228 for (j = 0; j < lenx; j++) {
2229 a_elem[0] = a_i[aij];
2230 a_elem[1] = a_i[aij + 1];
2231
2232 x_elem[0] = x_i[jx];
2233 x_elem[1] = x_i[jx + 1];
2234 {
2235 double head_e1, tail_e1;
2236 double d1;
2237 double d2;
2238 /* Real part */
2239 d1 = (double) a_elem[0] * x_elem[0];
2240 d2 = (double) -a_elem[1] * x_elem[1];
2241 {
2242 /* Compute double-double = double + double. */
2243 double e, t1, t2;
2244
2245 /* Knuth trick. */
2246 t1 = d1 + d2;
2247 e = t1 - d1;
2248 t2 = ((d2 - e) + (d1 - (t1 - e)));
2249
2250 /* The result is t1 + t2, after normalization. */
2251 head_e1 = t1 + t2;
2252 tail_e1 = t2 - (head_e1 - t1);
2253 }
2254 head_prod[0] = head_e1;
2255 tail_prod[0] = tail_e1;
2256 /* imaginary part */
2257 d1 = (double) a_elem[0] * x_elem[1];
2258 d2 = (double) a_elem[1] * x_elem[0];
2259 {
2260 /* Compute double-double = double + double. */
2261 double e, t1, t2;
2262
2263 /* Knuth trick. */
2264 t1 = d1 + d2;
2265 e = t1 - d1;
2266 t2 = ((d2 - e) + (d1 - (t1 - e)));
2267
2268 /* The result is t1 + t2, after normalization. */
2269 head_e1 = t1 + t2;
2270 tail_e1 = t2 - (head_e1 - t1);
2271 }
2272 head_prod[1] = head_e1;
2273 tail_prod[1] = tail_e1;
2274 }
2275 {
2276 double head_t, tail_t;
2277 double head_a, tail_a;
2278 double head_b, tail_b;
2279 /* Real part */
2280 head_a = head_sum[0];
2281 tail_a = tail_sum[0];
2282 head_b = head_prod[0];
2283 tail_b = tail_prod[0];
2284 {
2285 /* Compute double-double = double-double + double-double. */
2286 double bv;
2287 double s1, s2, t1, t2;
2288
2289 /* Add two hi words. */
2290 s1 = head_a + head_b;
2291 bv = s1 - head_a;
2292 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2293
2294 /* Add two lo words. */
2295 t1 = tail_a + tail_b;
2296 bv = t1 - tail_a;
2297 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2298
2299 s2 += t1;
2300
2301 /* Renormalize (s1, s2) to (t1, s2) */
2302 t1 = s1 + s2;
2303 s2 = s2 - (t1 - s1);
2304
2305 t2 += s2;
2306
2307 /* Renormalize (t1, t2) */
2308 head_t = t1 + t2;
2309 tail_t = t2 - (head_t - t1);
2310 }
2311 head_sum[0] = head_t;
2312 tail_sum[0] = tail_t;
2313 /* Imaginary part */
2314 head_a = head_sum[1];
2315 tail_a = tail_sum[1];
2316 head_b = head_prod[1];
2317 tail_b = tail_prod[1];
2318 {
2319 /* Compute double-double = double-double + double-double. */
2320 double bv;
2321 double s1, s2, t1, t2;
2322
2323 /* Add two hi words. */
2324 s1 = head_a + head_b;
2325 bv = s1 - head_a;
2326 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2327
2328 /* Add two lo words. */
2329 t1 = tail_a + tail_b;
2330 bv = t1 - tail_a;
2331 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2332
2333 s2 += t1;
2334
2335 /* Renormalize (s1, s2) to (t1, s2) */
2336 t1 = s1 + s2;
2337 s2 = s2 - (t1 - s1);
2338
2339 t2 += s2;
2340
2341 /* Renormalize (t1, t2) */
2342 head_t = t1 + t2;
2343 tail_t = t2 - (head_t - t1);
2344 }
2345 head_sum[1] = head_t;
2346 tail_sum[1] = tail_t;
2347 }
2348 aij += incaij;
2349 jx += incx;
2350 }
2351 {
2352 double cd[2];
2353 cd[0] = (double) alpha_i[0];
2354 cd[1] = (double) alpha_i[1];
2355 {
2356 /* Compute complex-extra = complex-extra * complex-double. */
2357 double head_a0, tail_a0;
2358 double head_a1, tail_a1;
2359 double head_t1, tail_t1;
2360 double head_t2, tail_t2;
2361 head_a0 = head_sum[0];
2362 tail_a0 = tail_sum[0];
2363 head_a1 = head_sum[1];
2364 tail_a1 = tail_sum[1];
2365 /* real part */
2366 {
2367 /* Compute double-double = double-double * double. */
2368 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
2369
2370 con = head_a0 * split;
2371 a11 = con - head_a0;
2372 a11 = con - a11;
2373 a21 = head_a0 - a11;
2374 con = cd[0] * split;
2375 b1 = con - cd[0];
2376 b1 = con - b1;
2377 b2 = cd[0] - b1;
2378
2379 c11 = head_a0 * cd[0];
2380 c21 =
2381 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
2382
2383 c2 = tail_a0 * cd[0];
2384 t1 = c11 + c2;
2385 t2 = (c2 - (t1 - c11)) + c21;
2386
2387 head_t1 = t1 + t2;
2388 tail_t1 = t2 - (head_t1 - t1);
2389 }
2390 {
2391 /* Compute double-double = double-double * double. */
2392 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
2393
2394 con = head_a1 * split;
2395 a11 = con - head_a1;
2396 a11 = con - a11;
2397 a21 = head_a1 - a11;
2398 con = cd[1] * split;
2399 b1 = con - cd[1];
2400 b1 = con - b1;
2401 b2 = cd[1] - b1;
2402
2403 c11 = head_a1 * cd[1];
2404 c21 =
2405 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
2406
2407 c2 = tail_a1 * cd[1];
2408 t1 = c11 + c2;
2409 t2 = (c2 - (t1 - c11)) + c21;
2410
2411 head_t2 = t1 + t2;
2412 tail_t2 = t2 - (head_t2 - t1);
2413 }
2414 head_t2 = -head_t2;
2415 tail_t2 = -tail_t2;
2416 {
2417 /* Compute double-double = double-double + double-double. */
2418 double bv;
2419 double s1, s2, t1, t2;
2420
2421 /* Add two hi words. */
2422 s1 = head_t1 + head_t2;
2423 bv = s1 - head_t1;
2424 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
2425
2426 /* Add two lo words. */
2427 t1 = tail_t1 + tail_t2;
2428 bv = t1 - tail_t1;
2429 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
2430
2431 s2 += t1;
2432
2433 /* Renormalize (s1, s2) to (t1, s2) */
2434 t1 = s1 + s2;
2435 s2 = s2 - (t1 - s1);
2436
2437 t2 += s2;
2438
2439 /* Renormalize (t1, t2) */
2440 head_t1 = t1 + t2;
2441 tail_t1 = t2 - (head_t1 - t1);
2442 }
2443 head_tmp1[0] = head_t1;
2444 tail_tmp1[0] = tail_t1;
2445 /* imaginary part */
2446 {
2447 /* Compute double-double = double-double * double. */
2448 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
2449
2450 con = head_a1 * split;
2451 a11 = con - head_a1;
2452 a11 = con - a11;
2453 a21 = head_a1 - a11;
2454 con = cd[0] * split;
2455 b1 = con - cd[0];
2456 b1 = con - b1;
2457 b2 = cd[0] - b1;
2458
2459 c11 = head_a1 * cd[0];
2460 c21 =
2461 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
2462
2463 c2 = tail_a1 * cd[0];
2464 t1 = c11 + c2;
2465 t2 = (c2 - (t1 - c11)) + c21;
2466
2467 head_t1 = t1 + t2;
2468 tail_t1 = t2 - (head_t1 - t1);
2469 }
2470 {
2471 /* Compute double-double = double-double * double. */
2472 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
2473
2474 con = head_a0 * split;
2475 a11 = con - head_a0;
2476 a11 = con - a11;
2477 a21 = head_a0 - a11;
2478 con = cd[1] * split;
2479 b1 = con - cd[1];
2480 b1 = con - b1;
2481 b2 = cd[1] - b1;
2482
2483 c11 = head_a0 * cd[1];
2484 c21 =
2485 (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
2486
2487 c2 = tail_a0 * cd[1];
2488 t1 = c11 + c2;
2489 t2 = (c2 - (t1 - c11)) + c21;
2490
2491 head_t2 = t1 + t2;
2492 tail_t2 = t2 - (head_t2 - t1);
2493 }
2494 {
2495 /* Compute double-double = double-double + double-double. */
2496 double bv;
2497 double s1, s2, t1, t2;
2498
2499 /* Add two hi words. */
2500 s1 = head_t1 + head_t2;
2501 bv = s1 - head_t1;
2502 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
2503
2504 /* Add two lo words. */
2505 t1 = tail_t1 + tail_t2;
2506 bv = t1 - tail_t1;
2507 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
2508
2509 s2 += t1;
2510
2511 /* Renormalize (s1, s2) to (t1, s2) */
2512 t1 = s1 + s2;
2513 s2 = s2 - (t1 - s1);
2514
2515 t2 += s2;
2516
2517 /* Renormalize (t1, t2) */
2518 head_t1 = t1 + t2;
2519 tail_t1 = t2 - (head_t1 - t1);
2520 }
2521 head_tmp1[1] = head_t1;
2522 tail_tmp1[1] = tail_t1;
2523 }
2524
2525 }
2526 y_elem[0] = y_i[iy];
2527 y_elem[1] = y_i[iy + 1];
2528 {
2529 double head_e1, tail_e1;
2530 double d1;
2531 double d2;
2532 /* Real part */
2533 d1 = (double) y_elem[0] * beta_i[0];
2534 d2 = (double) -y_elem[1] * beta_i[1];
2535 {
2536 /* Compute double-double = double + double. */
2537 double e, t1, t2;
2538
2539 /* Knuth trick. */
2540 t1 = d1 + d2;
2541 e = t1 - d1;
2542 t2 = ((d2 - e) + (d1 - (t1 - e)));
2543
2544 /* The result is t1 + t2, after normalization. */
2545 head_e1 = t1 + t2;
2546 tail_e1 = t2 - (head_e1 - t1);
2547 }
2548 head_tmp2[0] = head_e1;
2549 tail_tmp2[0] = tail_e1;
2550 /* imaginary part */
2551 d1 = (double) y_elem[0] * beta_i[1];
2552 d2 = (double) y_elem[1] * beta_i[0];
2553 {
2554 /* Compute double-double = double + double. */
2555 double e, t1, t2;
2556
2557 /* Knuth trick. */
2558 t1 = d1 + d2;
2559 e = t1 - d1;
2560 t2 = ((d2 - e) + (d1 - (t1 - e)));
2561
2562 /* The result is t1 + t2, after normalization. */
2563 head_e1 = t1 + t2;
2564 tail_e1 = t2 - (head_e1 - t1);
2565 }
2566 head_tmp2[1] = head_e1;
2567 tail_tmp2[1] = tail_e1;
2568 }
2569 {
2570 double head_t, tail_t;
2571 double head_a, tail_a;
2572 double head_b, tail_b;
2573 /* Real part */
2574 head_a = head_tmp1[0];
2575 tail_a = tail_tmp1[0];
2576 head_b = head_tmp2[0];
2577 tail_b = tail_tmp2[0];
2578 {
2579 /* Compute double-double = double-double + double-double. */
2580 double bv;
2581 double s1, s2, t1, t2;
2582
2583 /* Add two hi words. */
2584 s1 = head_a + head_b;
2585 bv = s1 - head_a;
2586 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2587
2588 /* Add two lo words. */
2589 t1 = tail_a + tail_b;
2590 bv = t1 - tail_a;
2591 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2592
2593 s2 += t1;
2594
2595 /* Renormalize (s1, s2) to (t1, s2) */
2596 t1 = s1 + s2;
2597 s2 = s2 - (t1 - s1);
2598
2599 t2 += s2;
2600
2601 /* Renormalize (t1, t2) */
2602 head_t = t1 + t2;
2603 tail_t = t2 - (head_t - t1);
2604 }
2605 head_tmp1[0] = head_t;
2606 tail_tmp1[0] = tail_t;
2607 /* Imaginary part */
2608 head_a = head_tmp1[1];
2609 tail_a = tail_tmp1[1];
2610 head_b = head_tmp2[1];
2611 tail_b = tail_tmp2[1];
2612 {
2613 /* Compute double-double = double-double + double-double. */
2614 double bv;
2615 double s1, s2, t1, t2;
2616
2617 /* Add two hi words. */
2618 s1 = head_a + head_b;
2619 bv = s1 - head_a;
2620 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2621
2622 /* Add two lo words. */
2623 t1 = tail_a + tail_b;
2624 bv = t1 - tail_a;
2625 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2626
2627 s2 += t1;
2628
2629 /* Renormalize (s1, s2) to (t1, s2) */
2630 t1 = s1 + s2;
2631 s2 = s2 - (t1 - s1);
2632
2633 t2 += s2;
2634
2635 /* Renormalize (t1, t2) */
2636 head_t = t1 + t2;
2637 tail_t = t2 - (head_t - t1);
2638 }
2639 head_tmp1[1] = head_t;
2640 tail_tmp1[1] = tail_t;
2641 }
2642 y_i[iy] = head_tmp1[0];
2643 y_i[iy + 1] = head_tmp1[1];
2644 ai += incai;
2645 iy += incy;
2646 }
2647 }
2648
2649 }
2650 }
2651
2652 FPU_FIX_STOP;
2653 }
2654 break;
2655 }
2656 }
2657