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