use milliseconds for timing

mpi_stencil2d_gt.cc

@@ -112,8 +112,7 @@ void set_rank_device(int n_ranks, int rank)
   gt::backend::clib::device_set(device);
 }
 
-// exchange in non-contiguous second dimension, staging into contiguous buffers
+// exchange in first dimension, staging into contiguous buffers on device
-// on device
 void boundary_exchange_x(MPI_Comm comm, int world_size, int rank,
                          gt::gtensor_device<double, 2>& d_z, int n_bnd,
                          bool stage_host = false)
@@ -362,8 +361,8 @@ int main(int argc, char** argv)
     d_dzdx_numeric = stencil2d_1d_5(d_z, stencil5) * scale;
     gt::synchronize();
   }
-  printf("%d/%d exchange time %0.8f\n", world_rank, world_size,
+  printf("%d/%d exchange time %0.8f ms\n", world_rank, world_size,
-         total_time / n_iter);
+         total_time / n_iter * 1000);
 
   gt::copy(d_dzdx_numeric, h_dzdx_numeric);
 

mpi_stencil2d_sycl.cc

@@ -14,6 +14,7 @@
  * by default.
  */
 
+#include <cassert>
 #include <cmath>
 #include <mpi.h>
 #include <stdio.h>
@@ -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,23 +49,24 @@ 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 out_idx = idx2(out_nrows, row, col);
-      int stride = nrows;
+      int in_base_idx = idx2(in_nrows, row, col);
-      out2d[base_idx] = (stencil5[0] * in2d[base_idx + 0 * stride] +
+      out2d[out_idx] = (stencil5[0] * in2d[in_base_idx + 0] +
-                         stencil5[1] * in2d[base_idx + 1 * stride] +
+                        stencil5[1] * in2d[in_base_idx + 1] +
-                         stencil5[2] * in2d[base_idx + 2 * stride] +
+                        stencil5[2] * in2d[in_base_idx + 2] +
-                         stencil5[3] * in2d[base_idx + 3 * stride] +
+                        stencil5[3] * in2d[in_base_idx + 3] +
-                         stencil5[4] * in2d[base_idx + 4 * stride]) *
+                        stencil5[4] * in2d[in_base_idx + 4]) *
                        scale;
     });
   });
@@ -62,41 +74,41 @@ auto stencil2d_1d_5(sycl::queue& q, int nrows, int ncols, double* out2d,
 }
 
 /*
- * Copy slice of second (non-contiguous) dimension of in array into contiguous
+ * Copy slice of first dimension of in array into contiguous
- * buffer out. In has dimension nrows x ncols, buf has dimension nrows x (end
+ * buffer out. In has dimension nrows x ncols, buf has dimension (end - start) x
- * -start + 1).
+ * ncols.
  */
-auto buf_from_view(sycl::queue& q, int nrows, double* buf, double* in,
+auto buf_from_view(sycl::queue& q, int ncols, int buf_nrows, double* buf,
-                   int start, int end)
+                   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
+ * Copy contiguous buffer into first dimension of array as a slice. Out has
- * slice. Out has dimension nrows x ncols, buf has dimension nrows x (end -
+ * dimension nrows x ncols, buf has dimension (end - start) * ncols.
- * start + 1). ncols is not actually needed or used, just need start to end
- * to be valid indexes.
  */
-auto buf_to_view(sycl::queue& q, int nrows, double* out, double* buf, int start,
+auto buf_to_view(sycl::queue& q, int ncols, int out_nrows, double* out,
-                 int end)
+                 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<double>(n * n, q);
+  int n_bnd = 2;
-  double* buf = sycl::malloc_host<double>(n * 2, q);
+  int n_with_ghost = n + 2 * n_bnd;
-  double* buf2 = sycl::malloc_host<double>(n * 2, q);
+  double* data = sycl::malloc_host<double>(n_with_ghost * n, q);
+  double* buf = sycl::malloc_host<double>(n_bnd * n, q);
+  double* buf2 = sycl::malloc_host<double>(n_bnd * n, q);
 
-  for (int i = 0; i < n; i++) {
   for (int j = 0; j < n; j++) {
-      data[j + i * n] = i + j / 1000.0;
+    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[idx2(n_bnd, 0, j)] = 100.0 + j;
-    buf2[i + n] = 100.0 + i + 0.1;
+    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 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 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();
+  buf_from_view(q, n, n_bnd, buf, n_with_ghost, data, 0, n_bnd).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]);
+    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 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
+// exchange in first dimension, staging into contiguous buffers on device
-// on device
+void boundary_exchange_x(MPI_Comm comm, int world_size, int rank,
-void boundary_exchange_y(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 z_size = n_local_with_ghost * n_global;
-  int dzdy_size = n_global * n_local;
+  int dzdx_size = n_local * n_global;
 
   double* h_z = sycl::malloc_host<double>(z_size, q);
   double* d_z = sycl::malloc_device<double>(z_size, q);
 
-  double* h_dzdy_numeric = sycl::malloc_host<double>(dzdy_size, q);
+  double* h_dzdx_numeric = sycl::malloc_host<double>(dzdx_size, q);
-  double* h_dzdy_actual = sycl::malloc_host<double>(dzdy_size, q);
+  double* h_dzdx_actual = sycl::malloc_host<double>(dzdx_size, q);
-  double* d_dzdy_numeric = sycl::malloc_device<double>(dzdy_size, q);
+  double* d_dzdx_numeric = sycl::malloc_device<double>(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 = [](double x, double y) { return x * x * x + y * y; };
-  auto fn_dzdy = [](double x, double y) { return 2 * x; };
+  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);
+    for (int i = 0; i < n_local; i++) {
-      h_dzdy_actual[idx2(n_global, j, i)] = fn_dzdy(xtmp, ytmp);
+      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 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 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,
+  printf("%d/%d exchange time %0.8f ms\n", world_rank, world_size,
-         total_time / n_iter);
+         total_time / n_iter * 1000);
 
-  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,
+    printf("%d la %f\n%d ln %f\n", world_rank, h_dzdx_actual[idx], world_rank,
-           h_dzdy_numeric[idx]);
+           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,
+    printf("%d ra %f\n%d rn %f\n", world_rank, h_dzdx_actual[idx], world_rank,
-           h_dzdy_numeric[idx]);
+           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);