Directly specify the halo width

[Goon83]

Following some examples to learn the halo function in DASH. One question is about how to specify the width of halo layer without StencilSpec. I found most example codes use the below steps to declare a array with halo.

   using StencilSpec_t = dash::halo::StencilSpec<StencilP_t, 4>;
   using HaloWrapper_t = dash::halo::HaloMatrixWrapper<Array_t>;

    StencilSpec_t stencil_spec(StencilP_t(-1, 0), StencilP_t(1, 0), StencilP_t(0, -1), StencilP_t(0, 1));
    HaloWrapper_t halo_array(data_array, stencil_spec);

After studying the API of HaloWrapper_t, it seems that the stencil_spec is the only way to ask for a halo layer.
https://codedocs.xyz/dash-project/dash/a01194.html

Thanks.
Bin

[dhinf]

How should a direct halo width specification look like?
Lets say, you want a two element halo width in each direction, does it mean including diagonal elements or not etc. . The thing is, you need later on the stencil spec to access the halo elements.

[dhinf]

Possible is a stencil spec generator for typical stencil pattern.
E.g. a full stencil with a specified width.
Maybe this could be something you are looking for.

[Goon83]

Hi dhinf,

Lets say, you want a two element halo width in each direction, does it mean including diagonal elements or not etc. . The thing is, you need later on the stencil spec to access the halo elements.

--> Yes, it includes diagonal elements. The general idea is to separate the halo spec and stencil operator. For a case with their own stencil implementation (such as our ArrayUDF), the DASH can also works as a storage layer. Anyway, I guess it is just an API effort. I assume that current DASH infers the width from the stencil spec (but I might be wrong). So, it can just expose this to user directly. In our ArrayUDF on HDF5 file, it provides the below API for a 2D array:

     std::vector<int> halo(2);
     halo[0] = 2, halo[1] = 2;
     std::vector<int> chunk(2);
     chunk[0] = 4, chunk[1] = 4;
     Array(chunk, halo, ...)
     The Array has a halo layer (/w length 2) for each direction. 
     Each chunk has size as 6 by 6 (except boundary chunk).

Possible is a stencil spec generator for typical stencil pattern.
E.g. a full stencil with a specified width.
Maybe this could be something you are looking for.

--> This could be a possible solution but as said above, the stencil width can be exposed to user.

[dhinf]

Hi Goon83,

The halo width is based on the stencil spec. But i think it should be possible to also provide the width. I'm interested in your use case. The halo wrapper takes the local part of a dash::matrix or narray (naming is a little bit misleading). Then it adds the halo part, which is a separate chunk of memory. How does ArrayUDF handles the halo memory?

Like i said, i'm really interested how you want to use DASH with the halo wrapper.

[Goon83]

Hi dhinf,
When data is stored in memory, ArrayUDF views the halo layer as non-halo data too. It basically merges halo and non-halo points as contiguous and normal chunk. ArrayUDF only knows the halo when it performs I/O, i.e., reading data into its memory.

I came up a simple and naive example (below) to simulate the request of ArrayUDF. It creates a 2D array on 4 MPI processing, giving the total size of 8 by 8 and the tile size of 4 by 4. Each process writes a subset of the 2D array (note: without halo). Then, each process read a subset of the array with halo layer. This example code works so far. But, I think the read with halo might be better served by the array created with halo layer. Otherwise, each access to halo layer may trigger a communication and thus bad performance.

Any thoughts on this? Thanks.

#include <unistd.h>
#include <iostream>
#include <cstddef>
#include <iomanip>

#include <libdash.h>
//#include <mpi.h>

using namespace std;

using std::cout;

int main(int argc, char *argv[])
{
  dash::init(&argc, &argv);

  size_t team_size = dash::Team::All().size();
  dash::TeamSpec<2> teamspec;
  teamspec.balance_extents();

  dash::global_unit_t myid = dash::myid();
  size_t num_units = dash::Team::All().size();

  if (num_units != 4)
  {
    cout << "Please run with mpirun -n 4" << endl;
    exit(-1);
  }

  size_t tilesize_x = 4;
  size_t tilesize_y = 4;
  size_t rows = tilesize_x * (num_units / 2);
  size_t cols = tilesize_y * (num_units / 2);
  dash::Matrix<int, 2> matrix(
      dash::SizeSpec<2>(
          rows,
          cols),
      dash::DistributionSpec<2>(),
      dash::Team::All(),
      teamspec);
  size_t matrix_size = rows * cols;
  DASH_ASSERT(matrix_size == matrix.size());
  DASH_ASSERT(rows == matrix.extent(0));
  DASH_ASSERT(cols == matrix.extent(1));

  if (0 == myid)
  {
    cout << "Matrix size: " << rows
         << " x " << cols
         << " == " << matrix_size
         << endl;
  }

  std::vector<size_t> my_start(2);
  std::vector<size_t> my_end(2);

  switch (myid)
  {
  case 0:
    my_start[0] = 0;
    my_start[1] = 0;
    break;
  case 1:
    my_start[0] = 0;
    my_start[1] = tilesize_y;
    break;
  case 2:
    my_start[0] = tilesize_x;
    my_start[1] = 0;
    break;
  case 3:
    my_start[0] = tilesize_x;
    my_start[1] = tilesize_y;
    break;
  default:
    break;
  }

  my_end[0] = my_start[0] + tilesize_x - 1;
  my_end[1] = my_start[1] + tilesize_y - 1;

  if (!myid)
    cout << "Assigning matrix values" << endl;

  cout << myid << "'s start = (" << my_start[0] << ", " << my_start[1] << "), my_end = " << my_end[0] << ", " << my_end[1] << " )" << endl;
  for (size_t i = my_start[0]; i <= my_end[0]; i++)
  {
    for (size_t k = my_start[1]; k <= my_end[1]; k++)
    {
      matrix[i][k] = myid;
    }
  }

  // Units waiting for value initialization
  dash::Team::All().barrier();

  // Read and assert values in matrix
  for (size_t i = my_start[0]; i <= my_end[0]; i++)
  {
    for (size_t k = my_start[1]; k <= my_end[1]; k++)
    {
      int value = matrix[i][k];
      int expected = myid;
      DASH_ASSERT(expected == value);
    }
  }

  if (!myid)
  {
    cout << "Print matrix values at rank 0: " << endl;
    for (size_t i = 0; i < matrix.extent(0); i++)
    {
      for (size_t k = 0; k < matrix.extent(1); k++)
      {
        int value = matrix[i][k];
        cout << value;
        cout << "   ";
      }
      cout << endl;
    }
  }

  std::vector<size_t> my_start_halo(2);
  std::vector<size_t> my_end_halo(2);

  switch (myid)
  {
  case 0:
    my_start_halo[0] = 0;
    my_start_halo[1] = 0;
    my_end_halo[0] = my_start_halo[0] + tilesize_x;
    my_end_halo[1] = my_start_halo[1] + tilesize_y;
    break;
  case 1:
    my_start_halo[0] = 0;
    my_start_halo[1] = tilesize_y - 1;
    my_end_halo[0] = my_start_halo[0] + tilesize_x;
    my_end_halo[1] = my_start_halo[1] + tilesize_y;
    break;
  case 2:
    my_start_halo[0] = tilesize_x - 1;
    my_start_halo[1] = 0;
    my_end_halo[0] = my_start_halo[0] + tilesize_x;
    my_end_halo[1] = my_start_halo[1] + tilesize_y;
    break;
  case 3:
    my_start_halo[0] = tilesize_x - 1;
    my_start_halo[1] = tilesize_y - 1;
    my_end_halo[0] = my_start_halo[0] + tilesize_x;
    my_end_halo[1] = my_start_halo[1] + tilesize_y;
    break;
  default:
    break;
  }

  sleep(myid * 3);

  cout << "myid =" << myid << "'s output (" << my_start_halo[0] << "," << my_start_halo[1] << ")->(" << my_end_halo[0] << "," << my_end_halo[1] << "):" << endl;
  // Read and assert values in matrix
  for (size_t i = my_start_halo[0]; i <= my_end_halo[0]; i++)
  {
    for (size_t k = my_start_halo[1]; k <= my_end_halo[1]; k++)
    {
      int value = matrix[i][k];
      cout << value << "  ";
    }
    cout << endl;
  }

  dash::Team::All().barrier();

  dash::finalize();
}

Bests,
Bin

[dhinf]

But, I think the read with halo might be better served by the array created with halo layer. Otherwise, each access to halo layer may trigger a communication and thus bad performance.

Indeed, each access on non local elements triggers a blocked communication call. So the performance is really bad. The halo wrapper takes the local part of the distributed Array and the stencil spec. Then it calculates the memory size required for the halo elements and requests new memory. At the moment it is not possible to use the subscript operator to access the halo elements. That's why you need the stencil iterator. But maybe a subscript operator based access layer for the halos can be added in future.

The idea behind the stencil operator is to distinguish between inner and boundary calculation to asynchronously load the halo elements while calculating the inner part. Also, we wanted to use the NArray with halo and without halo support (therefor we decided to use separate memory for the halo elements).

best
Denis

[Goon83]

Hi dhinf,
Thanks for clarification. I am a little confused with how to use stencil iterator to access halo elements. Could you point me to some sample code or document? It's more helpful if you can help to modify above example code to with halo and stencil iterator.

Thanks.
Bin

[dhinf]

Hi Bin,

here is an example:
https://github.com/dash-project/dash/blob/development/dash/examples/ex.02.matrix.halo.heat_equation/main.cpp

Within the iteration loop you can see how to use the iterator. If you use the feat-halo branch the comment in line 145 (slow version) isn't true anymore. The performance is more competitive in this branch.

[devreal]

@dhinf Can we merge this branch into development?

[Goon83]

hi dhinf,
Could you help to explain the blow code blocks? I am studying the slow version but have some difficult to locate right document to understand them. Maybe, I don't get the background of the code right.

Following the example code, it is a 2D matrix. The Stencil is defined as
StencilSpecT stencil_spec( StencilT(-1, 0), StencilT(1, 0), StencilT( 0, -1), StencilT(0, 1));
I assume this line forces HaloMatrixWrapperT to have 1 layer of halo at each direction when it is used here:

  HaloMatrixWrapperT halomat(matrix, bound_spec, stencil_spec);
  HaloMatrixWrapperT halomat2(matrix2, bound_spec, stencil_spec);

Then, the following two lines create a iterator/operator based on the stencil_spec again.

 auto stencil_op = halomat.stencil_operator(stencil_spec);
  auto stencil_op2 = halomat2.stencil_operator(stencil_spec);

So far, I think the Stencil should be 2D on a 2D data, i.e., accessing point at left/right/up/down of each point and do some calculation.

My confusing part is the below code block.

1. The iterator (value_at(...) ) on both inner/boundary are 1D-index based. How is the 2D Stencil mapped onto the 1D value_at based one? Any layout changes, or duplication exits?

2. Why does the separated too loop works? Not sure how the boundary points between inner and boundary are processed. Say, when an inner.value_at(1) access a cell out of inner (within boundary), does DASH automatically deal it?

  // slow version
    auto it_end = current_op->inner.end();
    for(auto it = current_op->inner.begin(); it != it_end; ++it)
    {
      auto core = *it;
      auto dtheta = (it.value_at(0) + it.value_at(1) - 2 * core) / (dx * dx) +
                    (it.value_at(2) + it.value_at(3) - 2 * core) /(dy * dy);
      new_begin[it.lpos()] = core + k * dtheta * dt;
    }

    // Wait until all Halo updates ready
    current_halo->wait();

    // Calculation of boundary Halo elements
    auto it_bend = current_op->boundary.end();
    for (auto it = current_op->boundary.begin(); it != it_bend; ++it) {
      auto core = *it;
      double dtheta =
          (it.value_at(0) + it.value_at(1) - 2 * core) / (dx * dx) +
          (it.value_at(2) + it.value_at(3) - 2 * core) / (dy * dy);
      new_begin[it.lpos()] = core + k * dtheta * dt;
    }

[dhinf]

@dhinf Can we merge this branch into development?

I need to clean up some code. I think at the end of this week or next week we can merge the branch into development.

[dhinf]

1. The iterator (value_at(...) ) on both inner/boundary are 1D-index based. How is the 2D Stencil mapped onto the 1D value_at based one? Any layout changes, or duplication exits?

This is the stencil point position within the stencilspec.
StencilT(-1, 0) -> value_at(0)
StencilT(1, 0) -> value_at(1)
StencilT( 0, -1) -> value_at(2)
StencilT(0, 1) -> value_at(3)

you also can use -> value_at(StencilT(-1,0)) . Here the iterator tries to find the position of stencil point within the given stencil spec and then calls value_at(position). So you have an indirection and results in slightly less performance.
I recommend the value_at with the position as argument.

2\. Why does the separated too loop works? Not sure how the boundary points between inner and boundary are processed. Say, when an inner.value_at(1) access a cell out of inner (within boundary), does DASH automatically deal it?

Depending on the stencil width the iterator for the inner part only iterates over elements where all stencil points accesses non halo elements only . The boundary iterator iterates over all center elements with at least one halo element access. The method value_at in this case tests whether the stencil point is located in the halo or the inner part.

[Goon83]

Hi dhinf,
Thanks for explanations. For the second question, I am still confused. Say the below example with the size of 5 by 5. "0" and "A" are inner elements and non 0 are halo. Will the inner iterator only works on the element "A" ? Which loop in above code will deal with "0"?

1 2 3 4 5
6 0 0 0 7
8 0 A 0 9
10 0 0 0 11
11 12 13 14 15

BTW,

1. I still want to raise the same question again for your consideration, could we separate the stencil and halo layer support in NArray? and also provide a generic API to access both inner boundary element.

2)As previously discussed, the stencil spec generator may help the first case. Any idea (sample code) how to do that in DASH? Some of our use case can have over 100s layers of halo. A stencil spec generator may be a temporary resolution to it.

Bests,
Bin

[dhinf]

Hi Bin,

for(auto it = current_op->inner.begin(); it != it_end; ++it)

works on A

for (auto it = current_op->boundary.begin(); it != it_bend; ++it) {

works on all 0 elements

I will look into your other questions next week.

Best
Denis

[Goon83]

Hi Denis,
Thanks for the answers. Could you elaborate more details about the iteration order (for both inner and boundary)? the cell layout order (e.g., row major)in memory?

I may rephrase my question with the below code again. As the code says, it accesses halo elements as normal (inner) element through [] operator or value_at. Having a halo layer can boost the performance. However, the code must be re-written to use HaloMatrixWrapper and boundary/inner iterator. I think it is fine to use the HaloMatrixWrapper but the boundary/inner iterator become too complexed. Could we just add a value_at operator to HaloMatrixWrapper? This operator access both boundary and inner elements with global coordinate. But the value_at decides to choose which boundary/inner iterator to use to access data. Definitely, users can explicitly control the update the boundary elements. Thanks.

Bests,
Bin

#include <unistd.h>
#include <iostream>
#include <cstddef>
#include <iomanip>

#include <libdash.h>
//#include <mpi.h>

using namespace std;

using std::cout;

int main(int argc, char *argv[])
{
  dash::init(&argc, &argv);

  size_t team_size = dash::Team::All().size();
  dash::TeamSpec<2> teamspec;
  teamspec.balance_extents();

  dash::global_unit_t myid = dash::myid();
  size_t num_units = dash::Team::All().size();

  if (num_units != 4)
  {
    cout << "Please run with mpirun -n 4" << endl;
    exit(-1);
  }

  size_t tilesize_x = 4;
  size_t tilesize_y = 4;
  size_t rows = tilesize_x * (num_units / 2);
  size_t cols = tilesize_y * (num_units / 2);
  dash::Matrix<int, 2> matrix(
      dash::SizeSpec<2>(
          rows,
          cols),
      dash::DistributionSpec<2>(),
      dash::Team::All(),
      teamspec);
  size_t matrix_size = rows * cols;
  DASH_ASSERT(matrix_size == matrix.size());
  DASH_ASSERT(rows == matrix.extent(0));
  DASH_ASSERT(cols == matrix.extent(1));

  if (0 == myid)
  {
    cout << "Matrix size: " << rows
         << " x " << cols
         << " == " << matrix_size
         << endl;
  }

  std::vector<size_t> my_start(2);
  std::vector<size_t> my_end(2);

  switch (myid)
  {
  case 0:
    my_start[0] = 0;
    my_start[1] = 0;
    break;
  case 1:
    my_start[0] = 0;
    my_start[1] = tilesize_y;
    break;
  case 2:
    my_start[0] = tilesize_x;
    my_start[1] = 0;
    break;
  case 3:
    my_start[0] = tilesize_x;
    my_start[1] = tilesize_y;
    break;
  default:
    break;
  }

  my_end[0] = my_start[0] + tilesize_x - 1;
  my_end[1] = my_start[1] + tilesize_y - 1;

  if (!myid)
    cout << "Assigning matrix values" << endl;

  cout << myid << "'s start = (" << my_start[0] << ", " << my_start[1] << "), my_end = " << my_end[0] << ", " << my_end[1] << " )" << endl;
  for (size_t i = my_start[0]; i <= my_end[0]; i++)
  {
    for (size_t k = my_start[1]; k <= my_end[1]; k++)
    {
      matrix[i][k] = myid;
    }
  }

  // Units waiting for value initialization
  dash::Team::All().barrier();

  // Read and assert values in matrix
  for (size_t i = my_start[0]; i <= my_end[0]; i++)
  {
    for (size_t k = my_start[1]; k <= my_end[1]; k++)
    {
      int value = matrix[i][k];
      int expected = myid;
      DASH_ASSERT(expected == value);
    }
  }

  if (!myid)
  {
    cout << "Print matrix values at rank 0: " << endl;
    for (size_t i = 0; i < matrix.extent(0); i++)
    {
      for (size_t k = 0; k < matrix.extent(1); k++)
      {
        int value = matrix[i][k];
        cout << value;
        cout << "   ";
      }
      cout << endl;
    }
  }

  std::vector<size_t> my_start_halo(2);
  std::vector<size_t> my_end_halo(2);

  switch (myid)
  {
  case 0:
    my_start_halo[0] = 0;
    my_start_halo[1] = 0;
    my_end_halo[0] = my_start_halo[0] + tilesize_x;
    my_end_halo[1] = my_start_halo[1] + tilesize_y;
    break;
  case 1:
    my_start_halo[0] = 0;
    my_start_halo[1] = tilesize_y - 1;
    my_end_halo[0] = my_start_halo[0] + tilesize_x;
    my_end_halo[1] = my_start_halo[1] + tilesize_y;
    break;
  case 2:
    my_start_halo[0] = tilesize_x - 1;
    my_start_halo[1] = 0;
    my_end_halo[0] = my_start_halo[0] + tilesize_x;
    my_end_halo[1] = my_start_halo[1] + tilesize_y;
    break;
  case 3:
    my_start_halo[0] = tilesize_x - 1;
    my_start_halo[1] = tilesize_y - 1;
    my_end_halo[0] = my_start_halo[0] + tilesize_x;
    my_end_halo[1] = my_start_halo[1] + tilesize_y;
    break;
  default:
    break;
  }

  sleep(myid * 3);

  cout << "myid =" << myid << "'s output (" << my_start_halo[0] << "," << my_start_halo[1] << ")->(" << my_end_halo[0] << "," << my_end_halo[1] << "):" << endl;
  // Read and assert values in matrix
  for (size_t i = my_start_halo[0]; i <= my_end_halo[0]; i++)
  {
    for (size_t k = my_start_halo[1]; k <= my_end_halo[1]; k++)
    {
      int value = matrix[i][k];
      cout << value << "  ";
    }
    cout << endl;
  }

  dash::Team::All().barrier();

  dash::finalize();
}