cg.c 10 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443
  1. /* StarPU --- Runtime system for heterogeneous multicore architectures.
  2. *
  3. * Copyright (C) 2010-2012 Université de Bordeaux 1
  4. *
  5. * StarPU is free software; you can redistribute it and/or modify
  6. * it under the terms of the GNU Lesser General Public License as published by
  7. * the Free Software Foundation; either version 2.1 of the License, or (at
  8. * your option) any later version.
  9. *
  10. * StarPU is distributed in the hope that it will be useful, but
  11. * WITHOUT ANY WARRANTY; without even the implied warranty of
  12. * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  13. *
  14. * See the GNU Lesser General Public License in COPYING.LGPL for more details.
  15. */
  16. #include <math.h>
  17. #include <assert.h>
  18. #include <sys/time.h>
  19. #include <starpu.h>
  20. #include <common/blas.h>
  21. #ifdef STARPU_USE_CUDA
  22. #include <cuda.h>
  23. #include <cublas.h>
  24. #endif
  25. #define FPRINTF(ofile, fmt, ...) do { if (!getenv("STARPU_SSILENT")) {fprintf(ofile, fmt, ## __VA_ARGS__); }} while(0)
  26. /*
  27. * Conjugate Gradient
  28. *
  29. * Input:
  30. * - matrix A
  31. * - vector b
  32. * - vector x (starting value)
  33. * - int i_max, error tolerance eps < 1.
  34. * Ouput:
  35. * - vector x
  36. *
  37. * Pseudo code:
  38. *
  39. * i <- 0
  40. * r <- b - Ax
  41. * d <- r
  42. * delta_new <- dot(r,r)
  43. * delta_0 <- delta_new
  44. *
  45. * while (i < i_max && delta_new > eps^2 delta_0)
  46. * {
  47. * q <- Ad
  48. * alpha <- delta_new/dot(d, q)
  49. * x <- x + alpha d
  50. *
  51. * If (i is divisible by 50)
  52. * r <- b - Ax
  53. * else
  54. * r <- r - alpha q
  55. *
  56. * delta_old <- delta_new
  57. * delta_new <- dot(r,r)
  58. * beta <- delta_new/delta_old
  59. * d <- r + beta d
  60. * i <- i + 1
  61. * }
  62. *
  63. */
  64. #include "cg.h"
  65. static int long long n = 1024;
  66. static int nblocks = 8;
  67. static int use_reduction = 1;
  68. static starpu_data_handle_t A_handle, b_handle, x_handle;
  69. static TYPE *A, *b, *x;
  70. static int i_max = 4000;
  71. static double eps = (10e-14);
  72. static starpu_data_handle_t r_handle, d_handle, q_handle;
  73. static TYPE *r, *d, *q;
  74. static starpu_data_handle_t dtq_handle, rtr_handle;
  75. static TYPE dtq, rtr;
  76. extern struct starpu_codelet accumulate_variable_cl;
  77. extern struct starpu_codelet accumulate_vector_cl;
  78. extern struct starpu_codelet bzero_variable_cl;
  79. extern struct starpu_codelet bzero_vector_cl;
  80. /*
  81. * Generate Input data
  82. */
  83. static void generate_random_problem(void)
  84. {
  85. int i, j;
  86. starpu_malloc((void **)&A, n*n*sizeof(TYPE));
  87. starpu_malloc((void **)&b, n*sizeof(TYPE));
  88. starpu_malloc((void **)&x, n*sizeof(TYPE));
  89. assert(A && b && x);
  90. for (j = 0; j < n; j++)
  91. {
  92. b[j] = (TYPE)1.0;
  93. x[j] = (TYPE)0.0;
  94. /* We take Hilbert matrix that is not well conditionned but definite positive: H(i,j) = 1/(1+i+j) */
  95. for (i = 0; i < n; i++)
  96. {
  97. A[n*j + i] = (TYPE)(1.0/(1.0+i+j));
  98. }
  99. }
  100. /* Internal vectors */
  101. starpu_malloc((void **)&r, n*sizeof(TYPE));
  102. starpu_malloc((void **)&d, n*sizeof(TYPE));
  103. starpu_malloc((void **)&q, n*sizeof(TYPE));
  104. assert(r && d && q);
  105. memset(r, 0, n*sizeof(TYPE));
  106. memset(d, 0, n*sizeof(TYPE));
  107. memset(q, 0, n*sizeof(TYPE));
  108. }
  109. static void free_data(void)
  110. {
  111. starpu_free(A);
  112. starpu_free(b);
  113. starpu_free(x);
  114. starpu_free(r);
  115. starpu_free(d);
  116. starpu_free(q);
  117. }
  118. static void register_data(void)
  119. {
  120. starpu_matrix_data_register(&A_handle, 0, (uintptr_t)A, n, n, n, sizeof(TYPE));
  121. starpu_vector_data_register(&b_handle, 0, (uintptr_t)b, n, sizeof(TYPE));
  122. starpu_vector_data_register(&x_handle, 0, (uintptr_t)x, n, sizeof(TYPE));
  123. starpu_vector_data_register(&r_handle, 0, (uintptr_t)r, n, sizeof(TYPE));
  124. starpu_vector_data_register(&d_handle, 0, (uintptr_t)d, n, sizeof(TYPE));
  125. starpu_vector_data_register(&q_handle, 0, (uintptr_t)q, n, sizeof(TYPE));
  126. starpu_variable_data_register(&dtq_handle, 0, (uintptr_t)&dtq, sizeof(TYPE));
  127. starpu_variable_data_register(&rtr_handle, 0, (uintptr_t)&rtr, sizeof(TYPE));
  128. if (use_reduction)
  129. {
  130. starpu_data_set_reduction_methods(q_handle, &accumulate_vector_cl, &bzero_vector_cl);
  131. starpu_data_set_reduction_methods(r_handle, &accumulate_vector_cl, &bzero_vector_cl);
  132. starpu_data_set_reduction_methods(dtq_handle, &accumulate_variable_cl, &bzero_variable_cl);
  133. starpu_data_set_reduction_methods(rtr_handle, &accumulate_variable_cl, &bzero_variable_cl);
  134. }
  135. }
  136. static void unregister_data(void)
  137. {
  138. starpu_data_unpartition(A_handle, 0);
  139. starpu_data_unpartition(b_handle, 0);
  140. starpu_data_unpartition(x_handle, 0);
  141. starpu_data_unpartition(r_handle, 0);
  142. starpu_data_unpartition(d_handle, 0);
  143. starpu_data_unpartition(q_handle, 0);
  144. starpu_data_unregister(A_handle);
  145. starpu_data_unregister(b_handle);
  146. starpu_data_unregister(x_handle);
  147. starpu_data_unregister(r_handle);
  148. starpu_data_unregister(d_handle);
  149. starpu_data_unregister(q_handle);
  150. starpu_data_unregister(dtq_handle);
  151. starpu_data_unregister(rtr_handle);
  152. }
  153. /*
  154. * Data partitioning filters
  155. */
  156. struct starpu_data_filter vector_filter;
  157. struct starpu_data_filter matrix_filter_1;
  158. struct starpu_data_filter matrix_filter_2;
  159. static void partition_data(void)
  160. {
  161. assert(n % nblocks == 0);
  162. /*
  163. * Partition the A matrix
  164. */
  165. /* Partition into contiguous parts */
  166. matrix_filter_1.filter_func = starpu_matrix_filter_block;
  167. matrix_filter_1.nchildren = nblocks;
  168. /* Partition into non-contiguous parts */
  169. matrix_filter_2.filter_func = starpu_matrix_filter_vertical_block;
  170. matrix_filter_2.nchildren = nblocks;
  171. /* A is in FORTRAN ordering, starpu_data_get_sub_data(A_handle, 2, i,
  172. * j) designates the block in column i and row j. */
  173. starpu_data_map_filters(A_handle, 2, &matrix_filter_1, &matrix_filter_2);
  174. /*
  175. * Partition the vectors
  176. */
  177. vector_filter.filter_func = starpu_vector_filter_block;
  178. vector_filter.nchildren = nblocks;
  179. starpu_data_partition(b_handle, &vector_filter);
  180. starpu_data_partition(x_handle, &vector_filter);
  181. starpu_data_partition(r_handle, &vector_filter);
  182. starpu_data_partition(d_handle, &vector_filter);
  183. starpu_data_partition(q_handle, &vector_filter);
  184. }
  185. /*
  186. * Debug
  187. */
  188. #if 0
  189. static void display_vector(starpu_data_handle_t handle, TYPE *ptr)
  190. {
  191. unsigned block_size = n / nblocks;
  192. unsigned b, ind;
  193. for (b = 0; b < nblocks; b++)
  194. {
  195. starpu_data_acquire(starpu_data_get_sub_data(handle, 1, b), STARPU_R);
  196. for (ind = 0; ind < block_size; ind++)
  197. {
  198. FPRINTF(stderr, "%2.2e ", ptr[b*block_size + ind]);
  199. }
  200. FPRINTF(stderr, "| ");
  201. starpu_data_release(starpu_data_get_sub_data(handle, 1, b));
  202. }
  203. FPRINTF(stderr, "\n");
  204. }
  205. static void display_matrix(void)
  206. {
  207. unsigned i, j;
  208. for (i = 0; i < n; i++)
  209. {
  210. for (j = 0; j < n; j++)
  211. {
  212. FPRINTF(stderr, "%2.2e ", A[j*n + i]);
  213. }
  214. FPRINTF(stderr, "\n");
  215. }
  216. }
  217. #endif
  218. /*
  219. * Main loop
  220. */
  221. static int cg(void)
  222. {
  223. double delta_new, delta_old, delta_0;
  224. double alpha, beta;
  225. int i = 0;
  226. int ret;
  227. /* r <- b */
  228. ret = copy_handle(r_handle, b_handle, nblocks);
  229. if (ret == -ENODEV) return ret;
  230. /* r <- r - A x */
  231. ret = gemv_kernel(r_handle, A_handle, x_handle, 1.0, -1.0, nblocks, use_reduction);
  232. if (ret == -ENODEV) return ret;
  233. /* d <- r */
  234. ret = copy_handle(d_handle, r_handle, nblocks);
  235. if (ret == -ENODEV) return ret;
  236. /* delta_new = dot(r,r) */
  237. ret = dot_kernel(r_handle, r_handle, rtr_handle, nblocks, use_reduction);
  238. if (ret == -ENODEV) return ret;
  239. starpu_data_acquire(rtr_handle, STARPU_R);
  240. delta_new = rtr;
  241. delta_0 = delta_new;
  242. starpu_data_release(rtr_handle);
  243. FPRINTF(stderr, "*************** INITIAL ************ \n");
  244. FPRINTF(stderr, "Delta 0: %e\n", delta_new);
  245. struct timeval start;
  246. struct timeval end;
  247. gettimeofday(&start, NULL);
  248. while ((i < i_max) && ((double)delta_new > (double)(eps*eps*delta_0)))
  249. {
  250. /* q <- A d */
  251. gemv_kernel(q_handle, A_handle, d_handle, 0.0, 1.0, nblocks, use_reduction);
  252. /* dtq <- dot(d,q) */
  253. dot_kernel(d_handle, q_handle, dtq_handle, nblocks, use_reduction);
  254. /* alpha = delta_new / dtq */
  255. starpu_data_acquire(dtq_handle, STARPU_R);
  256. alpha = delta_new/dtq;
  257. starpu_data_release(dtq_handle);
  258. /* x <- x + alpha d */
  259. axpy_kernel(x_handle, d_handle, alpha, nblocks);
  260. if ((i % 50) == 0)
  261. {
  262. /* r <- b */
  263. copy_handle(r_handle, b_handle, nblocks);
  264. /* r <- r - A x */
  265. gemv_kernel(r_handle, A_handle, x_handle, 1.0, -1.0, nblocks, use_reduction);
  266. }
  267. else
  268. {
  269. /* r <- r - alpha q */
  270. axpy_kernel(r_handle, q_handle, -alpha, nblocks);
  271. }
  272. /* delta_new = dot(r,r) */
  273. dot_kernel(r_handle, r_handle, rtr_handle, nblocks, use_reduction);
  274. starpu_data_acquire(rtr_handle, STARPU_R);
  275. delta_old = delta_new;
  276. delta_new = rtr;
  277. beta = delta_new / delta_old;
  278. starpu_data_release(rtr_handle);
  279. /* d <- beta d + r */
  280. scal_axpy_kernel(d_handle, beta, r_handle, 1.0, nblocks);
  281. if ((i % 10) == 0)
  282. {
  283. /* We here take the error as ||r||_2 / (n||b||_2) */
  284. double error = sqrt(delta_new/delta_0)/(1.0*n);
  285. FPRINTF(stderr, "*****************************************\n");
  286. FPRINTF(stderr, "iter %d DELTA %e - %e\n", i, delta_new, error);
  287. }
  288. i++;
  289. }
  290. gettimeofday(&end, NULL);
  291. double timing = (double)(((double)end.tv_sec - (double)start.tv_sec)*10e6 + ((double)end.tv_usec - (double)start.tv_usec));
  292. FPRINTF(stderr, "Total timing : %2.2f seconds\n", timing/10e6);
  293. FPRINTF(stderr, "Seconds per iteration : %2.2e\n", timing/10e6/i);
  294. return 0;
  295. }
  296. static int check(void)
  297. {
  298. return 0;
  299. }
  300. static void parse_args(int argc, char **argv)
  301. {
  302. int i;
  303. for (i = 1; i < argc; i++)
  304. {
  305. if (strcmp(argv[i], "-n") == 0)
  306. {
  307. n = (int long long)atoi(argv[++i]);
  308. continue;
  309. }
  310. if (strcmp(argv[i], "-maxiter") == 0)
  311. {
  312. i_max = atoi(argv[++i]);
  313. continue;
  314. }
  315. if (strcmp(argv[i], "-nblocks") == 0)
  316. {
  317. nblocks = atoi(argv[++i]);
  318. continue;
  319. }
  320. if (strcmp(argv[i], "-no-reduction") == 0)
  321. {
  322. use_reduction = 0;
  323. continue;
  324. }
  325. if (strcmp(argv[i], "-h") == 0)
  326. {
  327. FPRINTF(stderr, "usage: %s [-h] [-nblocks #blocks] [-n problem_size] [-no-reduction] [-maxiter i]\n", argv[0]);
  328. exit(-1);
  329. continue;
  330. }
  331. }
  332. }
  333. int main(int argc, char **argv)
  334. {
  335. int ret;
  336. #ifdef STARPU_QUICK_CHECK
  337. i_max = 16;
  338. #endif
  339. parse_args(argc, argv);
  340. ret = starpu_init(NULL);
  341. if (ret == -ENODEV)
  342. return 77;
  343. STARPU_CHECK_RETURN_VALUE(ret, "starpu_init");
  344. starpu_cublas_init();
  345. generate_random_problem();
  346. register_data();
  347. partition_data();
  348. ret = cg();
  349. if (ret == -ENODEV) goto enodev;
  350. ret = check();
  351. starpu_task_wait_for_all();
  352. unregister_data();
  353. free_data();
  354. starpu_cublas_shutdown();
  355. starpu_shutdown();
  356. return ret;
  357. enodev:
  358. starpu_shutdown();
  359. return 77;
  360. }