plu_implicit_example.c 8.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375
  1. /* StarPU --- Runtime system for heterogeneous multicore architectures.
  2. *
  3. * Copyright (C) 2010-2020 Université de Bordeaux, CNRS (LaBRI UMR 5800), Inria
  4. * Copyright (C) 2013 Thibaut Lambert
  5. *
  6. * StarPU is free software; you can redistribute it and/or modify
  7. * it under the terms of the GNU Lesser General Public License as published by
  8. * the Free Software Foundation; either version 2.1 of the License, or (at
  9. * your option) any later version.
  10. *
  11. * StarPU is distributed in the hope that it will be useful, but
  12. * WITHOUT ANY WARRANTY; without even the implied warranty of
  13. * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  14. *
  15. * See the GNU Lesser General Public License in COPYING.LGPL for more details.
  16. */
  17. #include <stdlib.h>
  18. #include <stdio.h>
  19. #include <string.h>
  20. #include <time.h>
  21. #include <math.h>
  22. #include <starpu.h>
  23. #include "pxlu.h"
  24. //#include "pxlu_kernels.h"
  25. #ifdef STARPU_HAVE_LIBNUMA
  26. #include <numaif.h>
  27. #endif
  28. static unsigned long size = 4096;
  29. static unsigned nblocks = 16;
  30. static unsigned check = 0;
  31. static int p = -1;
  32. static int q = -1;
  33. static unsigned display = 0;
  34. static unsigned no_prio = 0;
  35. #ifdef STARPU_HAVE_LIBNUMA
  36. static unsigned numa = 0;
  37. #endif
  38. static size_t allocated_memory = 0;
  39. static size_t allocated_memory_extra = 0;
  40. static starpu_data_handle_t *dataA_handles;
  41. static TYPE **dataA;
  42. int get_block_rank(unsigned i, unsigned j);
  43. static void parse_args(int argc, char **argv)
  44. {
  45. int i;
  46. for (i = 1; i < argc; i++)
  47. {
  48. if (strcmp(argv[i], "-size") == 0)
  49. {
  50. char *argptr;
  51. size = strtol(argv[++i], &argptr, 10);
  52. }
  53. if (strcmp(argv[i], "-nblocks") == 0)
  54. {
  55. char *argptr;
  56. nblocks = strtol(argv[++i], &argptr, 10);
  57. }
  58. if (strcmp(argv[i], "-check") == 0)
  59. {
  60. check = 1;
  61. }
  62. if (strcmp(argv[i], "-display") == 0)
  63. {
  64. display = 1;
  65. }
  66. if (strcmp(argv[i], "-numa") == 0)
  67. {
  68. #ifdef STARPU_HAVE_LIBNUMA
  69. numa = 1;
  70. #else
  71. fprintf(stderr, "Warning: libnuma is not available\n");
  72. #endif
  73. }
  74. if (strcmp(argv[i], "-p") == 0)
  75. {
  76. char *argptr;
  77. p = strtol(argv[++i], &argptr, 10);
  78. }
  79. if (strcmp(argv[i], "-q") == 0)
  80. {
  81. char *argptr;
  82. q = strtol(argv[++i], &argptr, 10);
  83. }
  84. if (strcmp(argv[i], "-h") == 0 || strcmp(argv[i], "-help") == 0 || strcmp(argv[i], "--help") == 0)
  85. {
  86. fprintf(stderr,"usage: %s [-size n] [-nblocks b] [-check] [-display] [-numa] [-p p] [-q q]\n", argv[0]);
  87. fprintf(stderr,"\np * q must be equal to the number of MPI nodes\n");
  88. exit(0);
  89. }
  90. }
  91. }
  92. unsigned STARPU_PLU(display_flag)(void)
  93. {
  94. return display;
  95. }
  96. static void fill_block_with_random(TYPE *blockptr, unsigned psize, unsigned pnblocks)
  97. {
  98. const unsigned block_size = (psize/pnblocks);
  99. unsigned i, j;
  100. for (i = 0; i < block_size; i++)
  101. for (j = 0; j < block_size; j++)
  102. {
  103. blockptr[j+i*block_size] = (TYPE)starpu_drand48();
  104. }
  105. }
  106. static void init_matrix(int rank)
  107. {
  108. #ifdef STARPU_HAVE_LIBNUMA
  109. if (numa)
  110. {
  111. fprintf(stderr, "Using INTERLEAVE policy\n");
  112. unsigned long nodemask = ((1<<0)|(1<<1));
  113. int ret = set_mempolicy(MPOL_INTERLEAVE, &nodemask, 3);
  114. if (ret)
  115. perror("set_mempolicy failed");
  116. }
  117. #endif
  118. /* Allocate a grid of data handles, not all of them have to be allocated later on */
  119. dataA_handles = calloc(nblocks*nblocks, sizeof(starpu_data_handle_t));
  120. dataA = calloc(nblocks*nblocks, sizeof(TYPE *));
  121. allocated_memory_extra += nblocks*nblocks*(sizeof(starpu_data_handle_t) + sizeof(TYPE *));
  122. size_t blocksize = (size_t)(size/nblocks)*(size/nblocks)*sizeof(TYPE);
  123. /* Allocate all the blocks that belong to this mpi node */
  124. unsigned long i,j;
  125. for (j = 0; j < nblocks; j++)
  126. {
  127. for (i = 0; i < nblocks; i++)
  128. {
  129. int block_rank = get_block_rank(i, j);
  130. TYPE **blockptr = &dataA[j+i*nblocks];
  131. // starpu_data_handle_t *handleptr = &dataA_handles[j+nblocks*i];
  132. starpu_data_handle_t *handleptr = &dataA_handles[j+nblocks*i];
  133. if (block_rank == rank)
  134. {
  135. /* This blocks should be treated by the current MPI process */
  136. /* Allocate and fill it */
  137. starpu_malloc((void **)blockptr, blocksize);
  138. allocated_memory += blocksize;
  139. //fprintf(stderr, "Rank %d : fill block (i = %d, j = %d)\n", rank, i, j);
  140. fill_block_with_random(*blockptr, size, nblocks);
  141. //fprintf(stderr, "Rank %d : fill block (i = %d, j = %d)\n", rank, i, j);
  142. if (i == j)
  143. {
  144. unsigned tmp;
  145. for (tmp = 0; tmp < size/nblocks; tmp++)
  146. {
  147. (*blockptr)[tmp*((size/nblocks)+1)] += (TYPE)10*nblocks;
  148. }
  149. }
  150. /* Register it to StarPU */
  151. starpu_matrix_data_register(handleptr, STARPU_MAIN_RAM,
  152. (uintptr_t)*blockptr, size/nblocks,
  153. size/nblocks, size/nblocks, sizeof(TYPE));
  154. }
  155. else
  156. {
  157. starpu_matrix_data_register(handleptr, -1,
  158. 0, size/nblocks,
  159. size/nblocks, size/nblocks, sizeof(TYPE));
  160. *blockptr = STARPU_POISON_PTR;
  161. }
  162. starpu_data_set_coordinates(*handleptr, 2, j, i);
  163. starpu_mpi_data_register(*handleptr, j+i*nblocks, block_rank);
  164. }
  165. }
  166. //display_all_blocks(nblocks, size/nblocks);
  167. }
  168. TYPE *STARPU_PLU(get_block)(unsigned i, unsigned j)
  169. {
  170. return dataA[j+i*nblocks];
  171. }
  172. int get_block_rank(unsigned i, unsigned j)
  173. {
  174. /* Take a 2D block cyclic distribution */
  175. /* NB: p (resp. q) is for "direction" i (resp. j) */
  176. return (j % q) * p + (i % p);
  177. }
  178. starpu_data_handle_t STARPU_PLU(get_block_handle)(unsigned i, unsigned j)
  179. {
  180. return dataA_handles[j+i*nblocks];
  181. }
  182. static void display_grid(int rank, unsigned pnblocks)
  183. {
  184. if (!display)
  185. return;
  186. //if (rank == 0)
  187. {
  188. fprintf(stderr, "2D grid layout (Rank %d): \n", rank);
  189. unsigned i, j;
  190. for (j = 0; j < pnblocks; j++)
  191. {
  192. for (i = 0; i < pnblocks; i++)
  193. {
  194. TYPE *blockptr = STARPU_PLU(get_block)(i, j);
  195. starpu_data_handle_t handle = STARPU_PLU(get_block_handle)(i, j);
  196. fprintf(stderr, "%d (data %p handle %p)", get_block_rank(i, j), blockptr, handle);
  197. }
  198. fprintf(stderr, "\n");
  199. }
  200. }
  201. }
  202. int main(int argc, char **argv)
  203. {
  204. int rank;
  205. int world_size;
  206. int ret;
  207. starpu_srand48((long int)time(NULL));
  208. parse_args(argc, argv);
  209. ret = starpu_mpi_init_conf(&argc, &argv, 1, MPI_COMM_WORLD, NULL);
  210. STARPU_CHECK_RETURN_VALUE(ret, "starpu_mpi_init_conf");
  211. starpu_mpi_comm_rank(MPI_COMM_WORLD, &rank);
  212. starpu_mpi_comm_size(MPI_COMM_WORLD, &world_size);
  213. if (p == -1 && q==-1)
  214. {
  215. fprintf(stderr, "Setting default values for p and q\n");
  216. p = (q % 2 == 0) ? 2 : 1;
  217. q = world_size / p;
  218. }
  219. STARPU_ASSERT_MSG(p*q == world_size, "p=%d, q=%d, world_size=%d\n", p, q, world_size);
  220. starpu_cublas_init();
  221. /*
  222. * Problem Init
  223. */
  224. init_matrix(rank);
  225. fprintf(stderr, "Rank %d: allocated (%d + %d) MB = %d MB\n", rank,
  226. (int)(allocated_memory/(1024*1024)),
  227. (int)(allocated_memory_extra/(1024*1024)),
  228. (int)((allocated_memory+allocated_memory_extra)/(1024*1024)));
  229. display_grid(rank, nblocks);
  230. TYPE *a_r = NULL;
  231. // STARPU_PLU(display_data_content)(a_r, size);
  232. if (check)
  233. {
  234. TYPE *x, *y;
  235. x = calloc(size, sizeof(TYPE));
  236. STARPU_ASSERT(x);
  237. y = calloc(size, sizeof(TYPE));
  238. STARPU_ASSERT(y);
  239. if (rank == 0)
  240. {
  241. unsigned ind;
  242. for (ind = 0; ind < size; ind++)
  243. x[ind] = (TYPE)starpu_drand48();
  244. }
  245. a_r = STARPU_PLU(reconstruct_matrix)(size, nblocks);
  246. if (rank == 0)
  247. STARPU_PLU(display_data_content)(a_r, size);
  248. // STARPU_PLU(compute_ax)(size, x, y, nblocks, rank);
  249. free(x);
  250. free(y);
  251. }
  252. double timing = STARPU_PLU(plu_main)(nblocks, rank, world_size, no_prio);
  253. /*
  254. * Report performance
  255. */
  256. if (rank == 0)
  257. {
  258. fprintf(stderr, "Computation took: %f ms\n", timing/1000);
  259. unsigned n = size;
  260. double flop = (2.0f*n*n*n)/3.0f;
  261. fprintf(stderr, "Synthetic GFlops : %2.2f\n", (flop/timing/1000.0f));
  262. }
  263. /*
  264. * Test Result Correctness
  265. */
  266. if (check)
  267. {
  268. /*
  269. * Compute || A - LU ||
  270. */
  271. STARPU_PLU(compute_lu_matrix)(size, nblocks, a_r);
  272. #if 0
  273. /*
  274. * Compute || Ax - LUx ||
  275. */
  276. unsigned ind;
  277. y2 = calloc(size, sizeof(TYPE));
  278. STARPU_ASSERT(y);
  279. if (rank == 0)
  280. {
  281. for (ind = 0; ind < size; ind++)
  282. {
  283. y2[ind] = (TYPE)0.0;
  284. }
  285. }
  286. STARPU_PLU(compute_lux)(size, x, y2, nblocks, rank);
  287. /* Compute y2 = y2 - y */
  288. CPU_AXPY(size, -1.0, y, 1, y2, 1);
  289. TYPE err = CPU_ASUM(size, y2, 1);
  290. int max = CPU_IAMAX(size, y2, 1);
  291. fprintf(stderr, "(A - LU)X Avg error : %e\n", err/(size*size));
  292. fprintf(stderr, "(A - LU)X Max error : %e\n", y2[max]);
  293. #endif
  294. }
  295. /*
  296. * Termination
  297. */
  298. starpu_cublas_shutdown();
  299. starpu_mpi_shutdown();
  300. return 0;
  301. }