You can not select more than 25 topics Topics must start with a chinese character,a letter or number, can include dashes ('-') and can be up to 35 characters long.

laswp_ncopy_8.c 7.3 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296
  1. /*********************************************************************/
  2. /* Copyright 2009, 2010 The University of Texas at Austin. */
  3. /* All rights reserved. */
  4. /* */
  5. /* Redistribution and use in source and binary forms, with or */
  6. /* without modification, are permitted provided that the following */
  7. /* conditions are met: */
  8. /* */
  9. /* 1. Redistributions of source code must retain the above */
  10. /* copyright notice, this list of conditions and the following */
  11. /* disclaimer. */
  12. /* */
  13. /* 2. Redistributions in binary form must reproduce the above */
  14. /* copyright notice, this list of conditions and the following */
  15. /* disclaimer in the documentation and/or other materials */
  16. /* provided with the distribution. */
  17. /* */
  18. /* THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF TEXAS AT */
  19. /* AUSTIN ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, */
  20. /* INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF */
  21. /* MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE */
  22. /* DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY OF TEXAS AT */
  23. /* AUSTIN OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, */
  24. /* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES */
  25. /* (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE */
  26. /* GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR */
  27. /* BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF */
  28. /* LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT */
  29. /* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT */
  30. /* OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE */
  31. /* POSSIBILITY OF SUCH DAMAGE. */
  32. /* */
  33. /* The views and conclusions contained in the software and */
  34. /* documentation are those of the authors and should not be */
  35. /* interpreted as representing official policies, either expressed */
  36. /* or implied, of The University of Texas at Austin. */
  37. /*********************************************************************/
  38. #include <stdio.h>
  39. #include "common.h"
  40. #define PREFETCHSIZE 4
  41. int CNAME(BLASLONG n, BLASLONG k1, BLASLONG k2, FLOAT *a, BLASLONG lda, blasint *ipiv, FLOAT *buffer){
  42. BLASLONG i, j, ip;
  43. blasint *piv;
  44. FLOAT *dx1, *dy1;
  45. FLOAT *dx2, *dy2;
  46. FLOAT *dx3, *dy3;
  47. FLOAT *dx4, *dy4;
  48. FLOAT *dx5, *dy5;
  49. FLOAT *dx6, *dy6;
  50. FLOAT *dx7, *dy7;
  51. FLOAT *dx8, *dy8;
  52. FLOAT atemp1, btemp1;
  53. FLOAT atemp2, btemp2;
  54. FLOAT atemp3, btemp3;
  55. FLOAT atemp4, btemp4;
  56. FLOAT atemp5, btemp5;
  57. FLOAT atemp6, btemp6;
  58. FLOAT atemp7, btemp7;
  59. FLOAT atemp8, btemp8;
  60. a--;
  61. ipiv += k1 - 1;
  62. if (n <= 0) return 0;
  63. if (k1 > k2) return 0;
  64. j = (n >> 3);
  65. if (j > 0) {
  66. do {
  67. piv = ipiv;
  68. i = k1;
  69. do {
  70. ip = *piv;
  71. piv ++;
  72. dx1 = a + i;
  73. dy1 = a + ip;
  74. dx2 = a + i + lda * 1;
  75. dy2 = a + ip + lda * 1;
  76. dx3 = a + i + lda * 2;
  77. dy3 = a + ip + lda * 2;
  78. dx4 = a + i + lda * 3;
  79. dy4 = a + ip + lda * 3;
  80. dx5 = a + i + lda * 4;
  81. dy5 = a + ip + lda * 4;
  82. dx6 = a + i + lda * 5;
  83. dy6 = a + ip + lda * 5;
  84. dx7 = a + i + lda * 6;
  85. dy7 = a + ip + lda * 6;
  86. dx8 = a + i + lda * 7;
  87. dy8 = a + ip + lda * 7;
  88. #ifdef __GNUC__
  89. __builtin_prefetch(dx1 + PREFETCHSIZE, 0, 1);
  90. __builtin_prefetch(dx2 + PREFETCHSIZE, 0, 1);
  91. __builtin_prefetch(dx3 + PREFETCHSIZE, 0, 1);
  92. __builtin_prefetch(dx4 + PREFETCHSIZE, 0, 1);
  93. __builtin_prefetch(dx5 + PREFETCHSIZE, 0, 1);
  94. __builtin_prefetch(dx6 + PREFETCHSIZE, 0, 1);
  95. __builtin_prefetch(dx7 + PREFETCHSIZE, 0, 1);
  96. __builtin_prefetch(dx8 + PREFETCHSIZE, 0, 1);
  97. #endif
  98. atemp1 = *dx1;
  99. btemp1 = *dy1;
  100. atemp2 = *dx2;
  101. btemp2 = *dy2;
  102. atemp3 = *dx3;
  103. btemp3 = *dy3;
  104. atemp4 = *dx4;
  105. btemp4 = *dy4;
  106. atemp5 = *dx5;
  107. btemp5 = *dy5;
  108. atemp6 = *dx6;
  109. btemp6 = *dy6;
  110. atemp7 = *dx7;
  111. btemp7 = *dy7;
  112. atemp8 = *dx8;
  113. btemp8 = *dy8;
  114. if (ip != i) {
  115. *dy1 = atemp1;
  116. *dy2 = atemp2;
  117. *dy3 = atemp3;
  118. *dy4 = atemp4;
  119. *dy5 = atemp5;
  120. *dy6 = atemp6;
  121. *dy7 = atemp7;
  122. *dy8 = atemp8;
  123. *(buffer + 0) = btemp1;
  124. *(buffer + 1) = btemp2;
  125. *(buffer + 2) = btemp3;
  126. *(buffer + 3) = btemp4;
  127. *(buffer + 4) = btemp5;
  128. *(buffer + 5) = btemp6;
  129. *(buffer + 6) = btemp7;
  130. *(buffer + 7) = btemp8;
  131. } else {
  132. *(buffer + 0) = atemp1;
  133. *(buffer + 1) = atemp2;
  134. *(buffer + 2) = atemp3;
  135. *(buffer + 3) = atemp4;
  136. *(buffer + 4) = atemp5;
  137. *(buffer + 5) = atemp6;
  138. *(buffer + 6) = atemp7;
  139. *(buffer + 7) = atemp8;
  140. }
  141. buffer += 8;
  142. i++;
  143. } while (i <= k2);
  144. a += 8 * lda;
  145. j --;
  146. } while (j > 0);
  147. }
  148. if (n & 4) {
  149. piv = ipiv;
  150. ip = *piv;
  151. piv ++;
  152. dx1 = a + k1;
  153. dy1 = a + ip;
  154. dx2 = a + k1 + lda * 1;
  155. dy2 = a + ip + lda * 1;
  156. dx3 = a + k1 + lda * 2;
  157. dy3 = a + ip + lda * 2;
  158. dx4 = a + k1 + lda * 3;
  159. dy4 = a + ip + lda * 3;
  160. i = k1;
  161. do {
  162. atemp1 = *dx1;
  163. atemp2 = *dx2;
  164. atemp3 = *dx3;
  165. atemp4 = *dx4;
  166. btemp1 = *dy1;
  167. btemp2 = *dy2;
  168. btemp3 = *dy3;
  169. btemp4 = *dy4;
  170. if (ip != i) {
  171. *dy1 = atemp1;
  172. *dy2 = atemp2;
  173. *dy3 = atemp3;
  174. *dy4 = atemp4;
  175. *(buffer + 0) = btemp1;
  176. *(buffer + 1) = btemp2;
  177. *(buffer + 2) = btemp3;
  178. *(buffer + 3) = btemp4;
  179. } else {
  180. *(buffer + 0) = atemp1;
  181. *(buffer + 1) = atemp2;
  182. *(buffer + 2) = atemp3;
  183. *(buffer + 3) = atemp4;
  184. }
  185. ip = *piv;
  186. piv ++;
  187. i++;
  188. dx1 = a + i;
  189. dy1 = a + ip;
  190. dx2 = a + i + lda * 1;
  191. dy2 = a + ip + lda * 1;
  192. dx3 = a + i + lda * 2;
  193. dy3 = a + ip + lda * 2;
  194. dx4 = a + i + lda * 3;
  195. dy4 = a + ip + lda * 3;
  196. buffer += 4;
  197. } while (i <= k2);
  198. a += 4 * lda;
  199. }
  200. if (n & 2) {
  201. piv = ipiv;
  202. i = k1;
  203. do {
  204. ip = *piv;
  205. piv ++;
  206. dx1 = a + i;
  207. dy1 = a + ip;
  208. dx2 = a + i + lda;
  209. dy2 = a + ip + lda;
  210. atemp1 = *dx1;
  211. btemp1 = *dy1;
  212. atemp2 = *dx2;
  213. btemp2 = *dy2;
  214. if (ip != i) {
  215. *dy1 = atemp1;
  216. *dy2 = atemp2;
  217. *(buffer + 0) = btemp1;
  218. *(buffer + 1) = btemp2;
  219. } else {
  220. *(buffer + 0) = atemp1;
  221. *(buffer + 1) = atemp2;
  222. }
  223. buffer += 2;
  224. i++;
  225. } while (i <= k2);
  226. a += 2 * lda;
  227. }
  228. if (n & 1) {
  229. piv = ipiv;
  230. i = k1;
  231. do {
  232. ip = *piv;
  233. piv ++;
  234. dx1 = a + i;
  235. dy1 = a + ip;
  236. atemp1 = *dx1;
  237. btemp1 = *dy1;
  238. if (ip != i) {
  239. *dy1 = atemp1;
  240. *buffer = btemp1;
  241. } else {
  242. *buffer = atemp1;
  243. }
  244. buffer ++;
  245. i++;
  246. } while (i <= k2);
  247. // a += lda;
  248. }
  249. return 0;
  250. }