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