diff --git a/mpi_stencil2d_sycl.cc b/mpi_stencil2d_sycl.cc index 073b493..8b18298 100644 --- a/mpi_stencil2d_sycl.cc +++ b/mpi_stencil2d_sycl.cc @@ -14,6 +14,7 @@ * by default. */ +#include #include #include #include @@ -22,6 +23,16 @@ #include "sycl/sycl.hpp" +inline void check(const char* fname, int line, int mpi_rval) +{ + if (mpi_rval != MPI_SUCCESS) { + printf("%s:%d error %d\n", fname, line, mpi_rval); + exit(2); + } +} + +#define CHECK(x) check(__FILE__, __LINE__, (x)) + static constexpr double stencil5[] = {1.0 / 12.0, -2.0 / 3.0, 0.0, 2.0 / 3.0, -1.0 / 12.0}; @@ -38,65 +49,66 @@ constexpr int idx2(int n, int row, int col) * * Returns sycl event, async with respect to host. */ -auto stencil2d_1d_5(sycl::queue& q, int nrows, int ncols, double* out2d, +auto stencil2d_1d_5(sycl::queue& q, int out_nrows, int out_ncols, double* out2d, const double* in2d, double scale) { // Note: swap index order; SYCL is row-major oriented, and this example // is col-major - auto range = sycl::range<2>(ncols, nrows); + int in_nrows = out_nrows + 4; + auto range = sycl::range<2>(out_ncols, out_nrows); auto e = q.submit([&](sycl::handler& cgh) { cgh.parallel_for(range, [=](sycl::item<2> item) { int row = item.get_id(1); int col = item.get_id(0); - int base_idx = idx2(nrows, row, col); - int stride = nrows; - out2d[base_idx] = (stencil5[0] * in2d[base_idx + 0 * stride] + - stencil5[1] * in2d[base_idx + 1 * stride] + - stencil5[2] * in2d[base_idx + 2 * stride] + - stencil5[3] * in2d[base_idx + 3 * stride] + - stencil5[4] * in2d[base_idx + 4 * stride]) * - scale; + int out_idx = idx2(out_nrows, row, col); + int in_base_idx = idx2(in_nrows, row, col); + out2d[out_idx] = (stencil5[0] * in2d[in_base_idx + 0] + + stencil5[1] * in2d[in_base_idx + 1] + + stencil5[2] * in2d[in_base_idx + 2] + + stencil5[3] * in2d[in_base_idx + 3] + + stencil5[4] * in2d[in_base_idx + 4]) * + scale; }); }); return e; } /* - * Copy slice of second (non-contiguous) dimension of in array into contiguous - * buffer out. In has dimension nrows x ncols, buf has dimension nrows x (end - * -start + 1). + * Copy slice of first dimension of in array into contiguous + * buffer out. In has dimension nrows x ncols, buf has dimension (end - start) x + * ncols. */ -auto buf_from_view(sycl::queue& q, int nrows, double* buf, double* in, - int start, int end) +auto buf_from_view(sycl::queue& q, int ncols, int buf_nrows, double* buf, + int in_nrows, double* in, int start, int end) { + assert(buf_nrows >= end - start); // Note: reverse index order because SYCL is row-major - auto range = sycl::range<2>(end - start, nrows); + auto range = sycl::range<2>(ncols, end - start); auto e = q.submit([&](sycl::handler& cgh) { cgh.parallel_for(range, [=](sycl::item<2> item) { int row = item.get_id(1); int col = item.get_id(0); - buf[idx2(nrows, row, col)] = in[idx2(nrows, row, start + col)]; + buf[idx2(buf_nrows, row, col)] = in[idx2(in_nrows, start + row, col)]; }); }); return e; } /* - * Copy contiguous buffer into second (non-contiguous) dimension of array as a - * slice. Out has dimension nrows x ncols, buf has dimension nrows x (end - - * start + 1). ncols is not actually needed or used, just need start to end - * to be valid indexes. + * Copy contiguous buffer into first dimension of array as a slice. Out has + * dimension nrows x ncols, buf has dimension (end - start) * ncols. */ -auto buf_to_view(sycl::queue& q, int nrows, double* out, double* buf, int start, - int end) +auto buf_to_view(sycl::queue& q, int ncols, int out_nrows, double* out, + int buf_nrows, double* buf, int start, int end) { + assert(buf_nrows >= end - start); // Note: reverse index order because SYCL is row-major - auto range = sycl::range<2>(end - start, nrows); + auto range = sycl::range<2>(ncols, end - start); auto e = q.submit([&](sycl::handler& cgh) { cgh.parallel_for(range, [=](sycl::item<2> item) { int row = item.get_id(1); int col = item.get_id(0); - out[idx2(nrows, row, start + col)] = buf[idx2(nrows, row, col)]; + out[idx2(out_nrows, start + row, col)] = buf[idx2(buf_nrows, row, col)]; }); }); return e; @@ -104,41 +116,43 @@ auto buf_to_view(sycl::queue& q, int nrows, double* out, double* buf, int start, void test_buf_view(sycl::queue& q, const int n) { - double* data = sycl::malloc_host(n * n, q); - double* buf = sycl::malloc_host(n * 2, q); - double* buf2 = sycl::malloc_host(n * 2, q); - - for (int i = 0; i < n; i++) { - for (int j = 0; j < n; j++) { - data[j + i * n] = i + j / 1000.0; + int n_bnd = 2; + int n_with_ghost = n + 2 * n_bnd; + double* data = sycl::malloc_host(n_with_ghost * n, q); + double* buf = sycl::malloc_host(n_bnd * n, q); + double* buf2 = sycl::malloc_host(n_bnd * n, q); + + for (int j = 0; j < n; j++) { + for (int i = 0; i < n_with_ghost; i++) { + data[idx2(n_with_ghost, i, j)] = (i - n_bnd) + j / 1000.0; } - buf2[i] = 100.0 + i; - buf2[i + n] = 100.0 + i + 0.1; + buf2[idx2(n_bnd, 0, j)] = 100.0 + j; + buf2[idx2(n_bnd, 1, j)] = 100.0 + j + 0.1; } - for (int i = 0; i < n; i++) { - for (int j = 0; j < n; j++) { - printf("data[%d, %d] = %f\n", j, i, data[j + i * n]); + for (int j = 0; j < n; j++) { + for (int i = 0; i < n; i++) { + printf("data[%d, %d] = %f\n", i, j, data[idx2(n_with_ghost, i, j)]); } } - for (int i = 0; i < 2; i++) { - for (int j = 0; j < n; j++) { - printf("buf2[%d, %d] = %f\n", j, i, buf2[j + i * n]); + for (int j = 0; j < n; j++) { + for (int i = 0; i < n_bnd; i++) { + printf("buf2[%d, %d] = %f\n", i, j, buf2[idx2(n_bnd, i, j)]); } } - buf_from_view(q, n, buf, data, 0, 2).wait(); - for (int i = 0; i < 2; i++) { - for (int j = 0; j < n; j++) { - printf("buf[%d, %d] = %f\n", j, i, buf[j + i * n]); + buf_from_view(q, n, n_bnd, buf, n_with_ghost, data, 0, n_bnd).wait(); + for (int j = 0; j < n; j++) { + for (int i = 0; i < n_bnd; i++) { + printf("buf[%d, %d] = %f\n", i, j, buf[idx2(n_bnd, i, j)]); } } - buf_to_view(q, n, data, buf2, n - 2, n).wait(); + buf_to_view(q, n, n_with_ghost, data, n_bnd, buf2, n - n_bnd, n).wait(); - for (int i = 0; i < n; i++) { - for (int j = 0; j < n; j++) { - printf("data[%d, %d] = %f\n", j, i, data[j + i * n]); + for (int j = 0; j < n; j++) { + for (int i = 0; i < n; i++) { + printf("data[%d, %d] = %f\n", i, j, data[idx2(n_with_ghost, i, j)]); } } } @@ -194,13 +208,13 @@ sycl::queue get_rank_queue(int n_ranks, int rank) cl::sycl::property::queue::in_order()}; } -// exchange in non-contiguous second dimension, staging into contiguous buffers -// on device -void boundary_exchange_y(MPI_Comm comm, int world_size, int rank, +// exchange in first dimension, staging into contiguous buffers on device +void boundary_exchange_x(MPI_Comm comm, int world_size, int rank, sycl::queue& q, int n_global, int n_local, int n_bnd, double* d_z, bool stage_host = false) { - int buf_size = n_global * n_bnd; + int n_local_with_ghost = n_local + 2 * n_bnd; + int buf_size = n_bnd * n_global; static double* sbuf_l = nullptr; static double* sbuf_r = nullptr; static double* rbuf_l = nullptr; @@ -234,7 +248,8 @@ void boundary_exchange_y(MPI_Comm comm, int world_size, int rank, if (rank_l >= 0) { // printf("rank_l = %d\n", rank_l); fflush(nullptr); // sbuf_l = d_z.view(_all, _s(n_bnd, 2 * n_bnd)); - auto e = buf_from_view(q, n_global, sbuf_l, d_z, n_bnd, 2 * n_bnd); + auto e = buf_from_view(q, n_global, n_bnd, sbuf_l, n_local_with_ghost, d_z, + n_bnd, 2 * n_bnd); if (stage_host) { q.copy(sbuf_l, h_sbuf_l, buf_size, e); /* @@ -251,7 +266,8 @@ void boundary_exchange_y(MPI_Comm comm, int world_size, int rank, if (rank_r < world_size) { // printf("rank_r = %d\n", rank_r); fflush(nullptr); // sbuf_r = d_z.view(_all, _s(-2 * n_bnd, -n_bnd)); - auto e = buf_from_view(q, n_global, sbuf_r, d_z, n_local, n_local + n_bnd); + auto e = buf_from_view(q, n_global, n_bnd, sbuf_r, n_local_with_ghost, d_z, + n_local, n_local + n_bnd); if (stage_host) { q.copy(sbuf_r, h_sbuf_r, buf_size, e); /* @@ -331,7 +347,7 @@ void boundary_exchange_y(MPI_Comm comm, int world_size, int rank, q.copy(h_rbuf_l, rbuf_l, buf_size); } // d_z.view(_all, _s(0, n_bnd)) = rbuf_l; - buf_to_view(q, n_global, d_z, rbuf_l, 0, n_bnd); + buf_to_view(q, n_global, n_local_with_ghost, d_z, n_bnd, rbuf_l, 0, n_bnd); } if (rank_r < world_size) { mpi_rval = MPI_Waitall(2, req_r, MPI_STATUSES_IGNORE); @@ -351,7 +367,8 @@ void boundary_exchange_y(MPI_Comm comm, int world_size, int rank, q.copy(h_rbuf_r, rbuf_r, buf_size); } // d_z.view(_all, _s(-n_bnd, _)) = rbuf_r; - buf_to_view(q, n_global, d_z, rbuf_r, n_local + n_bnd, n_local + 2 * n_bnd); + buf_to_view(q, n_global, n_local_with_ghost, d_z, n_bnd, rbuf_r, + n_local + n_bnd, n_local + 2 * n_bnd); } q.wait(); @@ -412,53 +429,53 @@ int main(int argc, char** argv) printf("stage_host = %d\n", stage_host); } - int z_size = n_global * n_local_with_ghost; - int dzdy_size = n_global * n_local; + int z_size = n_local_with_ghost * n_global; + int dzdx_size = n_local * n_global; double* h_z = sycl::malloc_host(z_size, q); double* d_z = sycl::malloc_device(z_size, q); - double* h_dzdy_numeric = sycl::malloc_host(dzdy_size, q); - double* h_dzdy_actual = sycl::malloc_host(dzdy_size, q); - double* d_dzdy_numeric = sycl::malloc_device(dzdy_size, q); + double* h_dzdx_numeric = sycl::malloc_host(dzdx_size, q); + double* h_dzdx_actual = sycl::malloc_host(dzdx_size, q); + double* d_dzdx_numeric = sycl::malloc_device(dzdx_size, q); double lx = 8; double dx = lx / n_global; double lx_local = lx / world_size; double scale = n_global / lx; - auto fn = [](double x, double y) { return x * x + y * y; }; - auto fn_dzdy = [](double x, double y) { return 2 * x; }; + auto fn = [](double x, double y) { return x * x * x + y * y; }; + auto fn_dzdx = [](double x, double y) { return 3 * x * x; }; struct timespec start, end; double iter_time = 0.0; double total_time = 0.0; double x_start = world_rank * lx_local; - for (int i = 0; i < n_local; i++) { - double xtmp = x_start + i * dx; - for (int j = 0; j < n_global; j++) { - double ytmp = j * dx; - h_z[idx2(n_global, j, i + n_bnd)] = fn(xtmp, ytmp); - h_dzdy_actual[idx2(n_global, j, i)] = fn_dzdy(xtmp, ytmp); + for (int j = 0; j < n_global; j++) { + double ytmp = j * dx; + for (int i = 0; i < n_local; i++) { + double xtmp = x_start + i * dx; + h_z[idx2(n_local_with_ghost, i + n_bnd, j)] = fn(xtmp, ytmp); + h_dzdx_actual[idx2(n_local, i, j)] = fn_dzdx(xtmp, ytmp); } } // fill boundary points on ends if (world_rank == 0) { - for (int i = 0; i < n_bnd; i++) { - double xtmp = (i - n_bnd) * dx; - for (int j = 0; j < n_global; j++) { - double ytmp = j * dx; - h_z[idx2(n_global, j, i)] = fn(xtmp, ytmp); + for (int j = 0; j < n_global; j++) { + double ytmp = j * dx; + for (int i = 0; i < n_bnd; i++) { + double xtmp = (i - n_bnd) * dx; + h_z[idx2(n_local_with_ghost, i, j)] = fn(xtmp, ytmp); } } } if (world_rank == world_size - 1) { - for (int i = 0; i < n_bnd; i++) { - double xtmp = lx + i * dx; - for (int j = 0; j < n_global; j++) { - double ytmp = j * dx; - h_z[idx2(n_global, j, n_bnd + n_local + i)] = fn(xtmp, ytmp); + for (int j = 0; j < n_global; j++) { + double ytmp = j * dx; + for (int i = 0; i < n_bnd; i++) { + double xtmp = lx + i * dx; + h_z[idx2(n_local_with_ghost, n_bnd + n_local + i, j)] = fn(xtmp, ytmp); } } } @@ -478,7 +495,7 @@ int main(int argc, char** argv) for (int i = 0; i < n_warmup + n_iter; i++) { clock_gettime(CLOCK_MONOTONIC, &start); - boundary_exchange_y(MPI_COMM_WORLD, world_size, world_rank, q, n_global, + boundary_exchange_x(MPI_COMM_WORLD, world_size, world_rank, q, n_global, n_local, n_bnd, d_z, stage_host); clock_gettime(CLOCK_MONOTONIC, &end); iter_time = @@ -489,28 +506,28 @@ int main(int argc, char** argv) } // do some calculation, to try to more closely simulate what happens in GENE - auto e = stencil2d_1d_5(q, n_global, n_local, d_dzdy_numeric, d_z, scale); + auto e = stencil2d_1d_5(q, n_local, n_global, d_dzdx_numeric, d_z, scale); e.wait(); } printf("%d/%d exchange time %0.8f\n", world_rank, world_size, total_time / n_iter); - q.copy(d_dzdy_numeric, h_dzdy_numeric, dzdy_size).wait(); + q.copy(d_dzdx_numeric, h_dzdx_numeric, dzdx_size).wait(); /* for (int i = 0; i < 5; i++) { int idx = idx2(n_global, 8, i); - printf("%d la %f\n%d ln %f\n", world_rank, h_dzdy_actual[idx], world_rank, - h_dzdy_numeric[idx]); + printf("%d la %f\n%d ln %f\n", world_rank, h_dzdx_actual[idx], world_rank, + h_dzdx_numeric[idx]); } for (int i = 0; i < 5; i++) { int idx = idx2(n_global, 8, n_local - 1 - i); - printf("%d ra %f\n%d rn %f\n", world_rank, h_dzdy_actual[idx], world_rank, - h_dzdy_numeric[idx]); + printf("%d ra %f\n%d rn %f\n", world_rank, h_dzdx_actual[idx], world_rank, + h_dzdx_numeric[idx]); } */ - double err_norm = diff_norm(q, dzdy_size, h_dzdy_numeric, h_dzdy_actual); + double err_norm = diff_norm(q, dzdx_size, h_dzdx_numeric, h_dzdx_actual); printf("%d/%d [%d:0x%08x] err_norm = %.8f\n", world_rank, world_size, device_id, vendor_id, err_norm);