1 /*********************************************************************/
2 /* */
3 /* Optimized BLAS libraries */
4 /* By Kazushige Goto <kgoto@tacc.utexas.edu> */
5 /* */
6 /* Copyright (c) The University of Texas, 2009. All rights reserved. */
7 /* UNIVERSITY EXPRESSLY DISCLAIMS ANY AND ALL WARRANTIES CONCERNING */
8 /* THIS SOFTWARE AND DOCUMENTATION, INCLUDING ANY WARRANTIES OF */
9 /* MERCHANTABILITY, FITNESS FOR ANY PARTICULAR PURPOSE, */
10 /* NON-INFRINGEMENT AND WARRANTIES OF PERFORMANCE, AND ANY WARRANTY */
11 /* THAT MIGHT OTHERWISE ARISE FROM COURSE OF DEALING OR USAGE OF */
12 /* TRADE. NO WARRANTY IS EITHER EXPRESS OR IMPLIED WITH RESPECT TO */
13 /* THE USE OF THE SOFTWARE OR DOCUMENTATION. */
14 /* Under no circumstances shall University be liable for incidental, */
15 /* special, indirect, direct or consequential damages or loss of */
16 /* profits, interruption of business, or related expenses which may */
17 /* arise from use of Software or Documentation, including but not */
18 /* limited to those resulting from defects in Software and/or */
19 /* Documentation, or loss or inaccuracy of data of any kind. */
20 /*********************************************************************/
21
22 #include <stdio.h>
23 #include "common.h"
24
25 #ifndef MINUS
26 #define a2 (a1 + 1)
27 #define a4 (a3 + 1)
28 #define a6 (a5 + 1)
29 #define a8 (a7 + 1)
30 #define a10 (a9 + 1)
31 #define a12 (a11 + 1)
32 #define a14 (a13 + 1)
33 #define a16 (a15 + 1)
34 #else
35 #define a2 (a1 - 1)
36 #define a4 (a3 - 1)
37 #define a6 (a5 - 1)
38 #define a8 (a7 - 1)
39 #define a10 (a9 - 1)
40 #define a12 (a11 - 1)
41 #define a14 (a13 - 1)
42 #define a16 (a15 - 1)
43 #endif
44
CNAME(BLASLONG n,BLASLONG k1,BLASLONG k2,FLOAT dummy1,FLOAT * a,BLASLONG lda,FLOAT * dummy2,BLASLONG dumy3,blasint * ipiv,BLASLONG incx)45 int CNAME(BLASLONG n, BLASLONG k1, BLASLONG k2, FLOAT dummy1, FLOAT *a, BLASLONG lda,
46 FLOAT *dummy2, BLASLONG dumy3, blasint *ipiv, BLASLONG incx){
47
48 BLASLONG i, j, ip1, ip2;
49 blasint *piv;
50 FLOAT *a1, *a3, *a5, *a7;
51 FLOAT *a9, *a11, *a13, *a15;
52 FLOAT *b1, *b2, *b3, *b4;
53 FLOAT *b5, *b6, *b7, *b8;
54 FLOAT *b9, *b10, *b11, *b12;
55 FLOAT *b13, *b14, *b15, *b16;
56 FLOAT A1, A2, B1, B2, A3, A4, B3, B4;
57 FLOAT A5, A6, B5, B6, A7, A8, B7, B8;
58 FLOAT A9, A10, B9, B10, A11, A12, B11, B12;
59 FLOAT A13, A14, B13, B14, A15, A16, B15, B16;
60
61 a--;
62 k1 --;
63
64 #ifndef MINUS
65 ipiv += k1;
66 #else
67 ipiv -= (k2 - 1) * incx;
68 #endif
69
70 if (n <= 0) return 0;
71
72 j = (n >> 3);
73 if (j > 0) {
74 do {
75 piv = ipiv;
76
77 #ifndef MINUS
78 a1 = a + k1 + 1;
79 #else
80 a1 = a + k2;
81 #endif
82
83 a3 = a1 + 1 * lda;
84 a5 = a1 + 2 * lda;
85 a7 = a1 + 3 * lda;
86 a9 = a1 + 4 * lda;
87 a11 = a1 + 5 * lda;
88 a13 = a1 + 6 * lda;
89 a15 = a1 + 7 * lda;
90
91 ip1 = *piv;
92 piv += incx;
93 ip2 = *piv;
94 piv += incx;
95
96 b1 = a + ip1;
97 b2 = a + ip2;
98
99 b3 = b1 + 1 * lda;
100 b4 = b2 + 1 * lda;
101 b5 = b1 + 2 * lda;
102 b6 = b2 + 2 * lda;
103 b7 = b1 + 3 * lda;
104 b8 = b2 + 3 * lda;
105
106 b9 = b1 + 4 * lda;
107 b10 = b2 + 4 * lda;
108 b11 = b1 + 5 * lda;
109 b12 = b2 + 5 * lda;
110 b13 = b1 + 6 * lda;
111 b14 = b2 + 6 * lda;
112 b15 = b1 + 7 * lda;
113 b16 = b2 + 7 * lda;
114
115 i = ((k2 - k1) >> 1);
116
117 if (i > 0) {
118 do {
119 B1 = *b1;
120 B2 = *b2;
121 B3 = *b3;
122 B4 = *b4;
123 B5 = *b5;
124 B6 = *b6;
125 B7 = *b7;
126 B8 = *b8;
127
128 B9 = *b9;
129 B10 = *b10;
130 B11 = *b11;
131 B12 = *b12;
132 B13 = *b13;
133 B14 = *b14;
134 B15 = *b15;
135 B16 = *b16;
136
137 A1 = *a1;
138 A2 = *a2;
139 A3 = *a3;
140 A4 = *a4;
141 A5 = *a5;
142 A6 = *a6;
143 A7 = *a7;
144 A8 = *a8;
145
146 A9 = *a9;
147 A10 = *a10;
148 A11 = *a11;
149 A12 = *a12;
150 A13 = *a13;
151 A14 = *a14;
152 A15 = *a15;
153 A16 = *a16;
154
155 ip1 = *piv;
156 piv += incx;
157 ip2 = *piv;
158 piv += incx;
159
160 if (b1 == a1) {
161 if (b2 == a1) {
162 *a1 = A2;
163 *a2 = A1;
164 *a3 = A4;
165 *a4 = A3;
166 *a5 = A6;
167 *a6 = A5;
168 *a7 = A8;
169 *a8 = A7;
170 *a9 = A10;
171 *a10 = A9;
172 *a11 = A12;
173 *a12 = A11;
174 *a13 = A14;
175 *a14 = A13;
176 *a15 = A16;
177 *a16 = A15;
178 } else
179 if (b2 != a2) {
180 *a2 = B2;
181 *b2 = A2;
182 *a4 = B4;
183 *b4 = A4;
184 *a6 = B6;
185 *b6 = A6;
186 *a8 = B8;
187 *b8 = A8;
188
189 *a10 = B10;
190 *b10 = A10;
191 *a12 = B12;
192 *b12 = A12;
193 *a14 = B14;
194 *b14 = A14;
195 *a16 = B16;
196 *b16 = A16;
197 }
198 } else
199 if (b1 == a2) {
200 if (b2 != a1) {
201 if (b2 == a2) {
202 *a1 = A2;
203 *a2 = A1;
204 *a3 = A4;
205 *a4 = A3;
206 *a5 = A6;
207 *a6 = A5;
208 *a7 = A8;
209 *a8 = A7;
210
211 *a9 = A10;
212 *a10 = A9;
213 *a11 = A12;
214 *a12 = A11;
215 *a13 = A14;
216 *a14 = A13;
217 *a15 = A16;
218 *a16 = A15;
219 } else {
220 *a1 = A2;
221 *a2 = B2;
222 *b2 = A1;
223 *a3 = A4;
224 *a4 = B4;
225 *b4 = A3;
226 *a5 = A6;
227 *a6 = B6;
228 *b6 = A5;
229 *a7 = A8;
230 *a8 = B8;
231 *b8 = A7;
232
233 *a9 = A10;
234 *a10 = B10;
235 *b10 = A9;
236 *a11 = A12;
237 *a12 = B12;
238 *b12 = A11;
239 *a13 = A14;
240 *a14 = B14;
241 *b14 = A13;
242 *a15 = A16;
243 *a16 = B16;
244 *b16 = A15;
245 }
246 }
247 } else {
248 if (b2 == a1) {
249 *a1 = A2;
250 *a2 = B1;
251 *b1 = A1;
252 *a3 = A4;
253 *a4 = B3;
254 *b3 = A3;
255 *a5 = A6;
256 *a6 = B5;
257 *b5 = A5;
258 *a7 = A8;
259 *a8 = B7;
260 *b7 = A7;
261
262 *a9 = A10;
263 *a10 = B9;
264 *b9 = A9;
265 *a11 = A12;
266 *a12 = B11;
267 *b11 = A11;
268 *a13 = A14;
269 *a14 = B13;
270 *b13 = A13;
271 *a15 = A16;
272 *a16 = B15;
273 *b15 = A15;
274 } else
275 if (b2 == a2) {
276 *a1 = B1;
277 *b1 = A1;
278 *a3 = B3;
279 *b3 = A3;
280 *a5 = B5;
281 *b5 = A5;
282 *a7 = B7;
283 *b7 = A7;
284
285 *a9 = B9;
286 *b9 = A9;
287 *a11 = B11;
288 *b11 = A11;
289 *a13 = B13;
290 *b13 = A13;
291 *a15 = B15;
292 *b15 = A15;
293 } else
294 if (b2 == b1) {
295 *a1 = B1;
296 *a2 = A1;
297 *b1 = A2;
298 *a3 = B3;
299 *a4 = A3;
300 *b3 = A4;
301 *a5 = B5;
302 *a6 = A5;
303 *b5 = A6;
304 *a7 = B7;
305 *a8 = A7;
306 *b7 = A8;
307
308 *a9 = B9;
309 *a10 = A9;
310 *b9 = A10;
311 *a11 = B11;
312 *a12 = A11;
313 *b11 = A12;
314 *a13 = B13;
315 *a14 = A13;
316 *b13 = A14;
317 *a15 = B15;
318 *a16 = A15;
319 *b15 = A16;
320 } else {
321 *a1 = B1;
322 *a2 = B2;
323 *b1 = A1;
324 *b2 = A2;
325 *a3 = B3;
326 *a4 = B4;
327 *b3 = A3;
328 *b4 = A4;
329 *a5 = B5;
330 *a6 = B6;
331 *b5 = A5;
332 *b6 = A6;
333 *a7 = B7;
334 *a8 = B8;
335 *b7 = A7;
336 *b8 = A8;
337
338 *a9 = B9;
339 *a10 = B10;
340 *b9 = A9;
341 *b10 = A10;
342 *a11 = B11;
343 *a12 = B12;
344 *b11 = A11;
345 *b12 = A12;
346 *a13 = B13;
347 *a14 = B14;
348 *b13 = A13;
349 *b14 = A14;
350 *a15 = B15;
351 *a16 = B16;
352 *b15 = A15;
353 *b16 = A16;
354 }
355 }
356
357 b1 = a + ip1;
358 b2 = a + ip2;
359
360 b3 = b1 + 1 * lda;
361 b4 = b2 + 1 * lda;
362 b5 = b1 + 2 * lda;
363 b6 = b2 + 2 * lda;
364 b7 = b1 + 3 * lda;
365 b8 = b2 + 3 * lda;
366
367 b9 = b1 + 4 * lda;
368 b10 = b2 + 4 * lda;
369 b11 = b1 + 5 * lda;
370 b12 = b2 + 5 * lda;
371 b13 = b1 + 6 * lda;
372 b14 = b2 + 6 * lda;
373 b15 = b1 + 7 * lda;
374 b16 = b2 + 7 * lda;
375
376 #ifndef MINUS
377 a1 += 2;
378 a3 += 2;
379 a5 += 2;
380 a7 += 2;
381 a9 += 2;
382 a11 += 2;
383 a13 += 2;
384 a15 += 2;
385 #else
386 a1 -= 2;
387 a3 -= 2;
388 a5 -= 2;
389 a7 -= 2;
390 a9 -= 2;
391 a11 -= 2;
392 a13 -= 2;
393 a15 -= 2;
394 #endif
395 i --;
396 } while (i > 0);
397 }
398
399 i = ((k2 - k1) & 1);
400
401 if (i > 0) {
402 A1 = *a1;
403 B1 = *b1;
404 A3 = *a3;
405 B3 = *b3;
406 A5 = *a5;
407 B5 = *b5;
408 A7 = *a7;
409 B7 = *b7;
410
411 A9 = *a9;
412 B9 = *b9;
413 A11 = *a11;
414 B11 = *b11;
415 A13 = *a13;
416 B13 = *b13;
417 A15 = *a15;
418 B15 = *b15;
419
420 *a1 = B1;
421 *b1 = A1;
422 *a3 = B3;
423 *b3 = A3;
424 *a5 = B5;
425 *b5 = A5;
426 *a7 = B7;
427 *b7 = A7;
428
429 *a9 = B9;
430 *b9 = A9;
431 *a11 = B11;
432 *b11 = A11;
433 *a13 = B13;
434 *b13 = A13;
435 *a15 = B15;
436 *b15 = A15;
437 }
438
439 a += 8 * lda;
440
441 j --;
442 } while (j > 0);
443 }
444
445 if (n & 4) {
446 piv = ipiv;
447
448 #ifndef MINUS
449 a1 = a + k1 + 1;
450 #else
451 a1 = a + k2;
452 #endif
453
454 a3 = a1 + 1 * lda;
455 a5 = a1 + 2 * lda;
456 a7 = a1 + 3 * lda;
457
458 ip1 = *piv;
459 piv += incx;
460 ip2 = *piv;
461 piv += incx;
462
463 b1 = a + ip1;
464 b2 = a + ip2;
465
466 b3 = b1 + 1 * lda;
467 b4 = b2 + 1 * lda;
468 b5 = b1 + 2 * lda;
469 b6 = b2 + 2 * lda;
470 b7 = b1 + 3 * lda;
471 b8 = b2 + 3 * lda;
472
473 i = ((k2 - k1) >> 1);
474
475 if (i > 0) {
476 do {
477 A1 = *a1;
478 A2 = *a2;
479 A3 = *a3;
480 A4 = *a4;
481 A5 = *a5;
482 A6 = *a6;
483 A7 = *a7;
484 A8 = *a8;
485
486 B1 = *b1;
487 B2 = *b2;
488 B3 = *b3;
489 B4 = *b4;
490 B5 = *b5;
491 B6 = *b6;
492 B7 = *b7;
493 B8 = *b8;
494
495 ip1 = *piv;
496 piv += incx;
497 ip2 = *piv;
498 piv += incx;
499
500 if (b1 == a1) {
501 if (b2 == a1) {
502 *a1 = A2;
503 *a2 = A1;
504 *a3 = A4;
505 *a4 = A3;
506 *a5 = A6;
507 *a6 = A5;
508 *a7 = A8;
509 *a8 = A7;
510 } else
511 if (b2 != a2) {
512 *a2 = B2;
513 *b2 = A2;
514 *a4 = B4;
515 *b4 = A4;
516 *a6 = B6;
517 *b6 = A6;
518 *a8 = B8;
519 *b8 = A8;
520 }
521 } else
522 if (b1 == a2) {
523 if (b2 != a1) {
524 if (b2 == a2) {
525 *a1 = A2;
526 *a2 = A1;
527 *a3 = A4;
528 *a4 = A3;
529 *a5 = A6;
530 *a6 = A5;
531 *a7 = A8;
532 *a8 = A7;
533 } else {
534 *a1 = A2;
535 *a2 = B2;
536 *b2 = A1;
537 *a3 = A4;
538 *a4 = B4;
539 *b4 = A3;
540 *a5 = A6;
541 *a6 = B6;
542 *b6 = A5;
543 *a7 = A8;
544 *a8 = B8;
545 *b8 = A7;
546 }
547 }
548 } else {
549 if (b2 == a1) {
550 *a1 = A2;
551 *a2 = B1;
552 *b1 = A1;
553 *a3 = A4;
554 *a4 = B3;
555 *b3 = A3;
556 *a5 = A6;
557 *a6 = B5;
558 *b5 = A5;
559 *a7 = A8;
560 *a8 = B7;
561 *b7 = A7;
562 } else
563 if (b2 == a2) {
564 *a1 = B1;
565 *b1 = A1;
566 *a3 = B3;
567 *b3 = A3;
568 *a5 = B5;
569 *b5 = A5;
570 *a7 = B7;
571 *b7 = A7;
572 } else
573 if (b2 == b1) {
574 *a1 = B1;
575 *a2 = A1;
576 *b1 = A2;
577 *a3 = B3;
578 *a4 = A3;
579 *b3 = A4;
580 *a5 = B5;
581 *a6 = A5;
582 *b5 = A6;
583 *a7 = B7;
584 *a8 = A7;
585 *b7 = A8;
586 } else {
587 *a1 = B1;
588 *a2 = B2;
589 *b1 = A1;
590 *b2 = A2;
591 *a3 = B3;
592 *a4 = B4;
593 *b3 = A3;
594 *b4 = A4;
595 *a5 = B5;
596 *a6 = B6;
597 *b5 = A5;
598 *b6 = A6;
599 *a7 = B7;
600 *a8 = B8;
601 *b7 = A7;
602 *b8 = A8;
603 }
604 }
605
606 b1 = a + ip1;
607 b2 = a + ip2;
608
609 b3 = b1 + 1 * lda;
610 b4 = b2 + 1 * lda;
611 b5 = b1 + 2 * lda;
612 b6 = b2 + 2 * lda;
613 b7 = b1 + 3 * lda;
614 b8 = b2 + 3 * lda;
615
616 #ifndef MINUS
617 a1 += 2;
618 a3 += 2;
619 a5 += 2;
620 a7 += 2;
621 #else
622 a1 -= 2;
623 a3 -= 2;
624 a5 -= 2;
625 a7 -= 2;
626 #endif
627 i --;
628 } while (i > 0);
629 }
630
631 i = ((k2 - k1) & 1);
632
633 if (i > 0) {
634 A1 = *a1;
635 B1 = *b1;
636 A3 = *a3;
637 B3 = *b3;
638 A5 = *a5;
639 B5 = *b5;
640 A7 = *a7;
641 B7 = *b7;
642
643 *a1 = B1;
644 *b1 = A1;
645 *a3 = B3;
646 *b3 = A3;
647 *a5 = B5;
648 *b5 = A5;
649 *a7 = B7;
650 *b7 = A7;
651 }
652
653 a += 4 * lda;
654 }
655
656 if (n & 2) {
657 piv = ipiv;
658
659 #ifndef MINUS
660 a1 = a + k1 + 1;
661 #else
662 a1 = a + k2;
663 #endif
664
665 a3 = a1 + 1 * lda;
666
667 ip1 = *piv;
668 piv += incx;
669 ip2 = *piv;
670 piv += incx;
671
672 b1 = a + ip1;
673 b2 = a + ip2;
674
675 b3 = b1 + 1 * lda;
676 b4 = b2 + 1 * lda;
677
678 i = ((k2 - k1) >> 1);
679
680 if (i > 0) {
681 do {
682 A1 = *a1;
683 A2 = *a2;
684 A3 = *a3;
685 A4 = *a4;
686
687 B1 = *b1;
688 B2 = *b2;
689 B3 = *b3;
690 B4 = *b4;
691
692 ip1 = *piv;
693 piv += incx;
694 ip2 = *piv;
695 piv += incx;
696
697 if (b1 == a1) {
698 if (b2 == a1) {
699 *a1 = A2;
700 *a2 = A1;
701 *a3 = A4;
702 *a4 = A3;
703 } else
704 if (b2 != a2) {
705 *a2 = B2;
706 *b2 = A2;
707 *a4 = B4;
708 *b4 = A4;
709 }
710 } else
711 if (b1 == a2) {
712 if (b2 != a1) {
713 if (b2 == a2) {
714 *a1 = A2;
715 *a2 = A1;
716 *a3 = A4;
717 *a4 = A3;
718 } else {
719 *a1 = A2;
720 *a2 = B2;
721 *b2 = A1;
722 *a3 = A4;
723 *a4 = B4;
724 *b4 = A3;
725 }
726 }
727 } else {
728 if (b2 == a1) {
729 *a1 = A2;
730 *a2 = B1;
731 *b1 = A1;
732 *a3 = A4;
733 *a4 = B3;
734 *b3 = A3;
735 } else
736 if (b2 == a2) {
737 *a1 = B1;
738 *b1 = A1;
739 *a3 = B3;
740 *b3 = A3;
741 } else
742 if (b2 == b1) {
743 *a1 = B1;
744 *a2 = A1;
745 *b1 = A2;
746 *a3 = B3;
747 *a4 = A3;
748 *b3 = A4;
749 } else {
750 *a1 = B1;
751 *a2 = B2;
752 *b1 = A1;
753 *b2 = A2;
754 *a3 = B3;
755 *a4 = B4;
756 *b3 = A3;
757 *b4 = A4;
758 }
759 }
760
761 b1 = a + ip1;
762 b2 = a + ip2;
763
764 b3 = b1 + 1 * lda;
765 b4 = b2 + 1 * lda;
766
767 #ifndef MINUS
768 a1 += 2;
769 a3 += 2;
770 #else
771 a1 -= 2;
772 a3 -= 2;
773 #endif
774 i --;
775 } while (i > 0);
776 }
777
778 i = ((k2 - k1) & 1);
779
780 if (i > 0) {
781 A1 = *a1;
782 B1 = *b1;
783 A3 = *a3;
784 B3 = *b3;
785 *a1 = B1;
786 *b1 = A1;
787 *a3 = B3;
788 *b3 = A3;
789 }
790
791 a += 2 * lda;
792 }
793
794 if (n & 1) {
795 piv = ipiv;
796
797 #ifndef MINUS
798 a1 = a + k1 + 1;
799 #else
800 a1 = a + k2;
801 #endif
802
803 ip1 = *piv;
804 piv += incx;
805 ip2 = *piv;
806 piv += incx;
807
808 b1 = a + ip1;
809 b2 = a + ip2;
810
811 i = ((k2 - k1) >> 1);
812
813 if (i > 0) {
814 do {
815 A1 = *a1;
816 A2 = *a2;
817 B1 = *b1;
818 B2 = *b2;
819
820 ip1 = *piv;
821 piv += incx;
822 ip2 = *piv;
823 piv += incx;
824
825 if (b1 == a1) {
826 if (b2 == a1) {
827 *a1 = A2;
828 *a2 = A1;
829 } else
830 if (b2 != a2) {
831 *a2 = B2;
832 *b2 = A2;
833 }
834 } else
835 if (b1 == a2) {
836 if (b2 != a1) {
837 if (b2 == a2) {
838 *a1 = A2;
839 *a2 = A1;
840 } else {
841 *a1 = A2;
842 *a2 = B2;
843 *b2 = A1;
844 }
845 }
846 } else {
847 if (b2 == a1) {
848 *a1 = A2;
849 *a2 = B1;
850 *b1 = A1;
851 } else
852 if (b2 == a2) {
853 *a1 = B1;
854 *b1 = A1;
855 } else
856 if (b2 == b1) {
857 *a1 = B1;
858 *a2 = A1;
859 *b1 = A2;
860 } else {
861 *a1 = B1;
862 *a2 = B2;
863 *b1 = A1;
864 *b2 = A2;
865 }
866 }
867
868 b1 = a + ip1;
869 b2 = a + ip2;
870
871 #ifndef MINUS
872 a1 += 2;
873 #else
874 a1 -= 2;
875 #endif
876 i --;
877 } while (i > 0);
878 }
879
880 i = ((k2 - k1) & 1);
881
882 if (i > 0) {
883 A1 = *a1;
884 B1 = *b1;
885 *a1 = B1;
886 *b1 = A1;
887 }
888 }
889
890 return 0;
891 }
892
893