cholesky.c 6.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264
  1. /* StarPU --- Runtime system for heterogeneous multicore architectures.
  2. *
  3. * Copyright (C) 2009-2011 Université de Bordeaux 1
  4. * Copyright (C) 2010 Mehdi Juhoor <mjuhoor@gmail.com>
  5. * Copyright (C) 2010, 2011, 2012 Centre National de la Recherche Scientifique
  6. *
  7. * StarPU is free software; you can redistribute it and/or modify
  8. * it under the terms of the GNU Lesser General Public License as published by
  9. * the Free Software Foundation; either version 2.1 of the License, or (at
  10. * your option) any later version.
  11. *
  12. * StarPU is distributed in the hope that it will be useful, but
  13. * WITHOUT ANY WARRANTY; without even the implied warranty of
  14. * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  15. *
  16. * See the GNU Lesser General Public License in COPYING.LGPL for more details.
  17. */
  18. #include "cholesky.h"
  19. #include "cholesky_kernels.h"
  20. #define __heap __attribute__ ((heap_allocated))
  21. /*
  22. * code to bootstrap the factorization
  23. * and construct the DAG
  24. */
  25. static void dw_cholesky(unsigned nblocks, unsigned size, unsigned ld, float *matA[nblocks][nblocks])
  26. {
  27. struct timeval start;
  28. struct timeval end;
  29. int x, y;
  30. /* create all the DAG nodes */
  31. unsigned i,j,k;
  32. for(x = 0; x < nblocks ; x++) {
  33. for (y = 0; y < nblocks; y++) {
  34. #pragma starpu register matA[x][y] size/nblocks*size/nblocks
  35. }
  36. }
  37. gettimeofday(&start, NULL);
  38. for (k = 0; k < nblocks; k++)
  39. {
  40. #ifdef STARPU_DEVEL
  41. # warning deal with prio and models
  42. #endif
  43. // int prio = STARPU_DEFAULT_PRIO;
  44. // if (!noprio) prio = STARPU_MAX_PRIO;
  45. chol_codelet_update_u11(matA[k][k], size/nblocks, ld);
  46. for (j = k+1; j<nblocks; j++)
  47. {
  48. // prio = STARPU_DEFAULT_PRIO;
  49. // if (!noprio&& (j == k+1)) prio = STARPU_MAX_PRIO;
  50. chol_codelet_update_u21(matA[k][k], matA[k][j], ld, ld, size/nblocks, size/nblocks);
  51. for (i = k+1; i<nblocks; i++)
  52. {
  53. if (i <= j)
  54. {
  55. // prio = STARPU_DEFAULT_PRIO;
  56. // if (!noprio && (i == k + 1) && (j == k +1) ) prio = STARPU_MAX_PRIO;
  57. chol_codelet_update_u22(matA[k][i],
  58. matA[k][j],
  59. matA[i][j],
  60. size/nblocks, size/nblocks, size/nblocks, ld, ld, ld);
  61. }
  62. }
  63. }
  64. }
  65. #pragma starpu wait
  66. for(x = 0; x < nblocks ; x++) {
  67. for (y = 0; y < nblocks; y++) {
  68. #pragma starpu unregister matA[x][y]
  69. }
  70. }
  71. gettimeofday(&end, NULL);
  72. double timing = (double)((end.tv_sec - start.tv_sec)*1000000 + (end.tv_usec - start.tv_usec));
  73. fprintf(stderr, "Computation took (in ms)\n");
  74. fprintf(stdout, "%2.2f\n", timing/1000);
  75. double flop = (1.0f*size*size*size)/3.0f;
  76. fprintf(stderr, "Synthetic GFlops : %2.2f\n", (flop/timing/1000.0f));
  77. }
  78. int main(int argc, char **argv)
  79. {
  80. /* create a simple definite positive symetric matrix example
  81. *
  82. * Hilbert matrix : h(i,j) = 1/(i+j+1)
  83. * */
  84. parse_args(argc, argv);
  85. #ifdef STARPU_DEVEL
  86. # warning todo
  87. #endif
  88. // struct starpu_conf conf;
  89. // starpu_conf_init(&conf);
  90. // conf.sched_policy_name = "heft";
  91. // conf.calibrate = 1;
  92. #pragma starpu initialize
  93. starpu_helper_cublas_init();
  94. float *bmat[nblocks][nblocks] __heap;
  95. unsigned i,j,x,y;
  96. for(x=0 ; x<nblocks ; x++)
  97. {
  98. for(y=0 ; y<nblocks ; y++)
  99. {
  100. starpu_malloc((void **)&bmat[x][y], BLOCKSIZE*BLOCKSIZE*sizeof(float));
  101. for (i = 0; i < BLOCKSIZE; i++)
  102. {
  103. for (j = 0; j < BLOCKSIZE; j++)
  104. {
  105. bmat[x][y][j +i*BLOCKSIZE] = (1.0f/(1.0f+(i+(x*BLOCKSIZE)+j+(y*BLOCKSIZE)))) + ((i+(x*BLOCKSIZE) == j+(y*BLOCKSIZE))?1.0f*size:0.0f);
  106. //mat[j +i*size] = ((i == j)?1.0f*size:0.0f);
  107. }
  108. }
  109. }
  110. }
  111. if (display) {
  112. for(y=0 ; y<nblocks ; y++)
  113. {
  114. for(x=0 ; x<nblocks ; x++)
  115. {
  116. printf("Block %d,%d :\n", x, y);
  117. for (j = 0; j < BLOCKSIZE; j++)
  118. {
  119. for (i = 0; i < BLOCKSIZE; i++)
  120. {
  121. if (i <= j) {
  122. printf("%2.2f\t", bmat[y][x][j +i*BLOCKSIZE]);
  123. }
  124. else {
  125. printf(".\t");
  126. }
  127. }
  128. printf("\n");
  129. }
  130. }
  131. }
  132. }
  133. dw_cholesky(nblocks, size, size/nblocks, bmat);
  134. if (display) {
  135. printf("Results:\n");
  136. for(y=0 ; y<nblocks ; y++)
  137. {
  138. for(x=0 ; x<nblocks ; x++)
  139. {
  140. printf("Block %d,%d :\n", x, y);
  141. for (j = 0; j < BLOCKSIZE; j++)
  142. {
  143. for (i = 0; i < BLOCKSIZE; i++)
  144. {
  145. if (i <= j) {
  146. printf("%2.2f\t", bmat[y][x][j +i*BLOCKSIZE]);
  147. }
  148. else {
  149. printf(".\t");
  150. }
  151. }
  152. printf("\n");
  153. }
  154. }
  155. }
  156. }
  157. float rmat[size * size] __heap;
  158. for(x=0 ; x<nblocks ; x++) {
  159. for(y=0 ; y<nblocks ; y++) {
  160. for (i = 0; i < BLOCKSIZE; i++) {
  161. for (j = 0; j < BLOCKSIZE; j++) {
  162. rmat[j+(y*BLOCKSIZE)+(i+(x*BLOCKSIZE))*size] = bmat[x][y][j +i*BLOCKSIZE];
  163. }
  164. }
  165. }
  166. }
  167. fprintf(stderr, "compute explicit LLt ...\n");
  168. for (j = 0; j < size; j++)
  169. {
  170. for (i = 0; i < size; i++)
  171. {
  172. if (i > j) {
  173. rmat[j+i*size] = 0.0f; // debug
  174. }
  175. }
  176. }
  177. float test_mat[size * size] __heap;
  178. SSYRK("L", "N", size, size, 1.0f,
  179. rmat, size, 0.0f, test_mat, size);
  180. fprintf(stderr, "comparing results ...\n");
  181. if (display) {
  182. for (j = 0; j < size; j++)
  183. {
  184. for (i = 0; i < size; i++)
  185. {
  186. if (i <= j) {
  187. printf("%2.2f\t", test_mat[j +i*size]);
  188. }
  189. else {
  190. printf(".\t");
  191. }
  192. }
  193. printf("\n");
  194. }
  195. }
  196. int correctness = 1;
  197. for(x = 0; x < nblocks ; x++)
  198. {
  199. for (y = 0; y < nblocks; y++)
  200. {
  201. for (i = (size/nblocks)*x ; i < (size/nblocks)*x+(size/nblocks); i++)
  202. {
  203. for (j = (size/nblocks)*y ; j < (size/nblocks)*y+(size/nblocks); j++)
  204. {
  205. if (i <= j)
  206. {
  207. float orig = (1.0f/(1.0f+i+j)) + ((i == j)?1.0f*size:0.0f);
  208. float err = abs(test_mat[j +i*size] - orig);
  209. if (err > 0.00001) {
  210. fprintf(stderr, "Error[%d, %d] --> %2.2f != %2.2f (err %2.2f)\n", i, j, test_mat[j +i*size], orig, err);
  211. correctness = 0;
  212. break;
  213. }
  214. }
  215. }
  216. }
  217. }
  218. }
  219. for(x=0 ; x<nblocks ; x++)
  220. {
  221. for(y=0 ; y<nblocks ; y++)
  222. {
  223. starpu_free((void *)bmat[x][y]);
  224. }
  225. }
  226. starpu_helper_cublas_shutdown();
  227. #pragma starpu shutdown
  228. assert(correctness);
  229. return 0;
  230. }