Saturday, 16 August 2014

What was that NDRange thing?

Digging a little deeper and pulling apart the OpenCL parts of the Mandelbrot Example

(This blog post refers to version 1.6.0 of the Brown Deer OpenCL libraries on the Parallella. If you are using newer versions there could be significant differences.)

I thought I'd stop and go back over the Mandelbrot example and make sure I understood how the thing worked. 

There were two things that seemed odd. Firstly, there was only one fork command (or rather a forka command) and I thought there would have been 16 and secondly, what was that clndrange_init1D command doing?

Turns out that the two are intimately linked. The ndRange controls the number of calls to your kernel and the fork kicks off the process. If you use forka you can pass in additional arguments. 

The nd part is short for n-Dimensional so the clndrange_init1D seems to suggest that the space created in the malloc statement as a 1-Dimensional space. Therefore, clndrange_init2D and clndrange_init3D would somehow create 2D and 3D spaces - or that is what I thought. It turns out that things are not exactly as clever as they might first appear.

Let's start from the begining. OpenCL has the standard call:

clEnqueueNDRangeKernel(cl_command_queue queue,
cl_kernel kernel,
cl_uint work_dims
const size_t *global_work_offset,
const size_t *global_work_size
const size_t *local_work_size,
cl_uint num_events, 
const cl_event *wait_list, 
cl_event *event);


Where:
  • work_dims is the number of dimensions (1, 2 or 3)
  • global_work_offset is the global ID offsets in each dimension
  • global_work_size is the number of work-items in each dimension
  • local_work_size is the number of work-items in a work-group, in each dimension

The last three are pointers to arrays of integers with for example global_work_offset[0] referring to the 1st dimension, global_work_offset[1] referring to the 2nd dimension etc.
(Note: this is the "official" version of what these are supposed to do)


The Brown Deer stdcl library replaces this one call with a choice of:

clndrange_init1d(gtoff0, gtsz0, ltsz0) OR
clndrange_init2d(gtoff0, gtsz0, ltsz0, gtoff1, gtsz1, ltsz1) OR
clndrange_init3d(gtoff0, gtsz0, ltsz0, gtoff1, gtsz1, ltsz1, gtoff2, gtsz2, ltsz2)

followed by:

clfork(context, devnum, kernel, ndrange, flags ) OR, if you want to pass additional arguments:
clforka(context, devnum, kernel, ndrange, flags [, arg0, ..., argN ]) 

where 

gtoff0 is the global_work_offset for dimension zero
gtsz0 is the global_work_size for dimension zero and
ltsz0 is the local_work_size for dimension zero etc.

So clndrange_init?d is a convenient way of setting those variables when is then passed into the fork command via the clndrange_t* variable.

So, what does it actually do?


clndrange_init1D


I started with a 1 dimensional experiment. After much digging and experimentation here's the DIRT!

global_work_offset does NOTHING, Seriously! Don't believe me... check out the khronos documentation here: http://www.khronos.org/registry/cl/sdk/1.0/docs/man/xhtml/clEnqueueNDRangeKernel.html
global_work_size is the number of times your kernel will be called (one call in opencl talk is one work-item) and
local_work_size is the number returned to you in your kernel when you call get_local_size(0)

That's it! No fancy partitioning of the data space. Nothing. It took me quite a while to realise that it was that simple. 

There are some quirks to this call. The combination of global_work_size and local_work_size are a bit sensitive. I have not tried every combination available but you often get a segmentation fault due to changes in both the global_work_size and local_work_size. For global_work_size less than 16 it seems to have problems if it and the local_work_size are not equal. Greater than 16 you could get an error "Exceeded maximum thread block size" but apart from getting an error, it seems to still work sometimes. It seems to work if global_work_size is a multiple of 16. If things are not right you can encounter Segmentation faults and Invalid Instruction faults or you may find that your kernel does not seem to be called at all. (I have used 16 here because I'm using an Epiphany-16. I should be referring to the number of cores on your accellerator.)

Let's have a look what the mandelbrot example did.

First, it allocated some space in the standard accelerator context stdacc with:

cl_uchar* pixels = (cl_uchar*) clmalloc(stdacc, width * height * 3, 0);

(The *3 is because it creates an RGB bit map which needs 3 bytes per pixel.)

After grabbing the compiled kernel from the context with clsym(...), it calls:

clndrange_t ndr = clndrange_init1d(0, height, 16);

so it wants is kernel to be called height times. The 16 is actually ignored because the kernel never calls get_local_size(0).

Finally it calls:

clforka(stdacc, 0, krn, &ndr, CL_EVENT_WAIT, iterations, width, startx, starty, dx_over_width, dy_over_height, pixels);

where it passes in the width of each line as an argument along with the data structure, pixels. The kernel then generates one horizontal line of the final picture using the call get_global_id(0) to determine which line it is on. (A slightly cleaner way would have been to pass width as the third argument in the clndrange_init1d() and then to call get_local_size(0)).

So much for clndrange_init1D. The mandelbrot set is a good example of a problem where the calculation of each point is independent of the next. If you have such a problem then this is a simple model that would be sufficient.

Next I tried a 2 dimensional version.


clndrange_init2D


Let me first say that I found clndrange_init2D to be a little temperamental. It seems that some combinations of values result in the kernel not being called at all. What's more, once it is not working, regressing to the previous state did not result it in working again. It was very frustrating. Therefore, everything written below must be read in with the knowledge that I gave up without actually understanding what was going on.

The 2D call seems (see previous paragraph) to call the kernel gtsz0 multiplied by gtsz1 times. The local size values ltsz0 and ltsz1 are merely returned by get_local_size(0|1).

clndrange_init2d( gtoff0, gtsz0, ltsz0, gtoff1, gtsz1, ltsz1)

gtoff0: NULL! Nothing to do here
gtsz0: the number of calls for the first dimension 
ltsz0: the value returned by calling get_local_size(0)
gtoff1: NULL, as you might expect
gtsz1: the number of calls for the second dimension (one for each gtsz0 call)
ltsz1: the number returning by calling get_local_size(1)

I didn't go into clndrange_init3D but I'd lay money on it working in the same was as clndrange_init2D.


nD Example 


I took the mandelbrot example and removed all the fancy mathematics. I allocated two chunks of global memory and wrote a 1D kernel and a 2D kernel that just initialised the space to a given value.

The critical bits are:

int  bytesPerCore = 16; // how many bytes we want each core to process
int workItems = 32;     // the total number workItems (threads sent to the Epiphany)

wrkArea1D = (cl_uchar*) clmalloc(stdacc, workItems * bytesPerCore, 0);
wrkArea2D = (cl_uchar*) clmalloc(stdacc, workItems * bytesPerCore, 0);

clndrange_t ndr1D = clndrange_init1d(NULL, workItems, bytesPerCore); 
clndrange_t ndr2D = clndrange_init2d(NULL, workItems, bytesPerCore/4, NULL, workItems, bytesPerCore/4);

and then in a common function call:

clforka(stdacc, 0, krn, ndr, CL_EVENT_WAIT, 1, rows, cols, wrkArea);

where 
krn is either "k_init1D" or "k_init2D"
ndr is either ndr1D or ndr2D and
wrkArea is either wrkArea1D or wrkArea2D

The kernel k_init1D works in the same way as the mandelbrot example processing 16 bytes of data in each of the 32 calls using get_global_id(0) to figure out where it should be.

The k_init2D kernel breaks the data set of the same size into 4x4 "tiles" so that adjacent data could be processed at the same time. I thought that I could cast the global data into a 2 dimensional array but that didn't work so I had to do all of the offset arithmetic in code. While this is not difficult it did make the 2D kernel considerably longer and given that speed is of the essence I would suggest that the only reason to do this would be if the algorithm overall would work better in 2D than in 1D (or that you dig unnecessary complexity).

In the resultant data sets I include the value of get_global_id(0) to show which call processed which chunk of data. The 1D data has the global id in the first column and the result of get_local_size(0) in the second column. The first five lines are:

0       16      1       1       1       1       1       1       1       1       1       1       1       1       1       1
1       16      1       1       1       1       1       1       1       1       1       1       1       1       1       1
2       16      1       1       1       1       1       1       1       1       1       1       1       1       1       1
3       16      1       1       1       1       1       1       1       1       1       1       1       1       1       1
4       16      1       1       1       1       1       1       1       1       1       1       1       1       1       1

The 2D data set includes the result of get_global_id(1) which shows it jumping around between 28 and 31. This is the value returned by get_global_id(1) as at the last call to the kernel. The top left value is the value returned by get_global_id(0) and the next value is value returned by get_global_id(1). The results from the first four tiles (1 to 3) are:

0       28     1       1       1       31     1       1       2       31     1       1       3       29     1       1
1       1       1       1       1       1       1       1       1       1       1       1       1       1       1       1
1       1       1       1       1       1       1       1       1       1       1       1       1       1       1       1
1       1       1       1       1       1       1       1       1       1       1       1       1       1       1       1

The 2D call was also a lot more sensitive than the 1D. As I mentioned above some combinations of values in the 2D version work and some don't. I suspect that the total number of calls must be a multiple of the number of cores available. So, before writing a big kernel with the assumption of clndrange_init2D working, think about how you want to process your data and check that a simple version works first and 

The nD example code


The code for the example is a little long to list here so grab it from github. From you eclipse workspace directory execute:

git clone https://github.com/nickoppen/nD.git

load it into eclipse and push, pull, tweak and generally rip it up till you are an ndrange EXPERT. And, if I've gotten anything wrong, please let me know.


Up Next


I don't know about you but I'm getting a bit tired of waiting around for Eclipse. AndyC posted a procedure on getting Code::Blocks working on the Parallella (http://forums.parallella.org/viewtopic.php?f=13&t=1658). 

I've installed Code::Blocks and it does perform better than Eclipse. Andy's procedure is designed for the standard SDK and I got into a mess when I tried to compile an OpenCL program on it so I'm going to have to do some tweaking. If I make some progress I'll write a post about it.

Sunday, 3 August 2014

OpenCL on the Parallella using Eclipse

Writing OpenCL code on the Parallella board using Eclipse for people who just want to write code and not become Linux experts.

Finally, what I would call my first "execution" blog post. I got my board a couple of months ago and while it's been fun playing with it, achieving the ultimate goal - writing code - has been frustrating to say the least.

I've tried a number of different approaches only to realise (yet again) that keeping it simple while you don't know what you are doing is by far the best approach. 

Where I ended up last week was - "how about I try and get some existing code running in eclipse on the parallella". Hopefully, the code will be compilable and correct so the only think I'll have to get right are the settings. To this end I found a short program that generated a colourised mandelbrot set posted by dar on the Parallella site. Don't try and find this code because it has some pretty glaring errors. If you want to try this procedure use the code below.

This really is a beginners guide. If you are already familiar with Eclipse you won't need most of this. I also don't make any claim about how efficient the code is.


0. Before you begin


The prerequisites are:

  • the Parallella SDK
  • the Brown Deer Technology SDK
  • the environment variables PATH and LD_LIBRARY_PATH set correctly
  • the Eclipse development environment


The first three come with the Ubuntu 14.04 release and I'm assuming future releases will also include them. For the record, the environment variables on my system are:

LD_LIBRARY_PATH=/usr/local/browndeer/lib:/usr/local/lib:/opt/adapteva/esdk/tools/host/lib:


PATH=/usr/local/browndeer/bin:/opt/adapteva/esdk/tools/host/bin:/opt/adapteva/esdk/tools/e-gnu/bin:/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin:/usr/games:/usr/local/games

Installing eclipse is a matter of:

sudo apt-get install eclipse

Once you down load eclipse you need to C/C++ development addin CDT. To get this you use Help | Get New Software and pack your patience - it is SLOOOOW.

On tricky thing: for a program to gain access to the Epiphany co-processor it must be run by root. That means that you must log in as root (using sudo won't work) and therefore you need a root password. To reset the root password you run:

sudo passwd

which will give you the usual new password prompt followed by the are you really sure prompt. Then to run you program you use the command:

su

from your non-root login (linaro by default). This will then prompt you for the root password. Once you have logged on as root, check the environment variables PATH and LD_LIBRARY_PATH. They are needed at run time so root must have them set as described above.

1. Your first Eclipse Project


1.1 Get Yourself a New Project

To start with you need a new C/C++ project. The File | New | Project pops up the box to choose what type of project you want. I always choose C++ but for this particular project a standard C Project would be fine. I called my project mandelbrot (not very original I know).

1.2 Set Up your Compile Settings

Not a lot to do here but absolutely critical.

Your compile settings are accessible from the Project | Properties menu. For Properties to be active you need to have the project tab and your project within it selected.

Your tool chain should look like this:



Note I've got the C++ compiler and linker in there.

Your includes should look like this:



This will tell the compiler to find all of the Brown Deer Technology stuff.

Finally, your linker settings should look like this:


I'm pretty sure that this is all you need. I did fumble around a lot with other settings so if this does not work please let me know.

1.3 Some Code

For CL projects you need host code and Epiphany code. The host code is compiled using the gcc/g++ compilers/linker etc in the tool chain. The Epiphany code is compiled at run time by the Brown Deer JIT compiler.

I created a source file folder (src) under the mandelbrot project folder. For the embedded paths in the code to work you should do this as well. The host code I ended up with was:

// The modifications porting this code to OpenCL are
// Copyright (c) 2012 Brown Deer Technology, LLC.
//
// Mandelbrot.c
// Written by User:Evercat
//
// This draws the Mandelbrot set and spits out a .bmp file.
// Should be quite portable (endian issues have been taken
// care of, for example)
//
// Released under the GNU Free Documentation License
// or the GNU Public License, whichever you prefer:
// 9 February, 2004.

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdcl.h>
#include <errno.h>

#define OUTFILE "./mandelbrot.bmp"

#define WIDTH 1024
#define HEIGHT 768

#define CENTRE_X -0.5
#define CENTRE_Y 0
#define ZOOM 300

#define ITERATIONS 1024  // Higher is more detailed, but slower...

// Plotting functions and parameters...

#define bailoutr(n) (5*n  )
#define bailoutg(n) (20*n  )
#define bailoutb(n) 0

#define min(a,b) ((a<b)?a:b)

// Colours for the set itself...

#define IN_SET_R 0
#define IN_SET_G 0
#define IN_SET_B 0

void drawbmp(int width, int height, unsigned char* pixels, char * filename);

/////////////////////////////////// MAIN PROGRAM ///////////////////////////////////

int main (void)
{

   float startx; float endx;
   float starty; float endy;
   float dx; float dy;
   float dx_over_width,dy_over_height;

   char kern[] = "../src/mandel_kern.cl";
   void * openHandle;

   int iterations = ITERATIONS;
   int width = WIDTH;
   int height = HEIGHT;

   char strInfo[20];
   FILE * pFile;

   pFile = fopen(kern, "r");
   if (pFile == NULL)
   {
  printf("Opening the Kernel file: %s produced an error(%d). Make sure that the source code variable kern has a valid path to the cl code and that the code is readable.\n", kern, errno);
  exit(0);
   }
   else
  fclose(pFile); // only open the file to check that it is there and readable

   pFile = fopen("./debug", "w");

   fprintf(pFile, "About to malloc pixels\n");
   cl_uchar* pixels = (cl_uchar*) clmalloc(stdacc, width * height * 3, 0);

   startx = CENTRE_X - ((float) WIDTH / (ZOOM * 2));
   endx = CENTRE_X + ((float) WIDTH / (ZOOM * 2));

   starty = CENTRE_Y - ((float) HEIGHT / (ZOOM * 2));
   endy = CENTRE_Y + ((float) HEIGHT / (ZOOM * 2));

   fprintf(pFile, "Plotting from (%f, %f) to (%f, %f)\n", startx, starty, endx, endy);

   dx = endx - startx;
   dy = endy - starty;
   dx_over_width = dx / width;
   dy_over_height = dy / height;


   fprintf(pFile, "Opening kernel file:%s\n", kern);
   openHandle = clopen(stdacc, kern, CLLD_NOW);

   fprintf(pFile, "Getting the kernel with clsym\n");
   cl_kernel krn = clsym(stdacc, openHandle, "mandel_kern", CLLD_NOW);

   clGetKernelInfo(krn, CL_KERNEL_FUNCTION_NAME, sizeof(strInfo), strInfo, NULL);
   fprintf(pFile, "The kernel is called: %s\n", strInfo);

   fprintf(pFile, "Calling clndrange\n");
   clndrange_t ndr = clndrange_init1d(0, height, 16);

   fprintf(pFile, "Calling clforka\n");
   clforka(stdacc, 0, krn, &ndr, CL_EVENT_WAIT,
      iterations, width, startx, starty, dx_over_width, dy_over_height, pixels);

   fprintf(pFile, "Transferring memory contents from the Epiphany using clmsync\n");
   clmsync(stdacc, 0, pixels, CL_MEM_HOST|CL_EVENT_WAIT);

   fprintf(pFile, "Calling drawbmp\n");
   drawbmp(width, height, pixels, OUTFILE);

   fprintf(pFile, "Saved bitmap to %s. Done.\n", OUTFILE);
   clfree(pixels);
   fclose(pFile);

   return 0;
}


void drawbmp (int width, int height, unsigned char* pixels, char * filename) {

   unsigned int headers[13];
   FILE * outfile;
   int extrabytes;
   int paddedsize;
   int x; int y; int n;

   extrabytes = 4 - ((width * 3) % 4); // How many bytes of padding to add to
                                       // each horizontal line - the size of
                                       // which must be a multiple of 4 bytes.
   if (extrabytes == 4)
      extrabytes = 0;

   paddedsize = ((width * 3) + extrabytes) * height;

   // Headers...

   headers[0]  = paddedsize + 54;      // bfSize (whole file size)
   headers[1]  = 0;                    // bfReserved (both)
   headers[2]  = 54;                   // bfOffbits
   headers[3]  = 40;                   // biSize
   headers[4]  = width;  // biWidth
   headers[5]  = height; // biHeight
                                       // 6 will be written directly...
   headers[7]  = 0;                    // biCompression
   headers[8]  = paddedsize;           // biSizeImage
   headers[9]  = 0;                    // biXPelsPerMeter
   headers[10] = 0;                    // biYPelsPerMeter
   headers[11] = 0;                    // biClrUsed
   headers[12] = 0;                    // biClrImportant

   outfile = fopen (filename, "wb");

   // Headers begin...
   // When printing ints and shorts, write out 1 character at time to
   // avoid endian issues.

   fprintf (outfile, "BM");

   for (n = 0; n <= 5; n++)
   {
      fprintf(outfile, "%c", headers[n] & 0x000000FF);
      fprintf(outfile, "%c", (headers[n] & 0x0000FF00) >> 8);
      fprintf(outfile, "%c", (headers[n] & 0x00FF0000) >> 16);
      fprintf(outfile, "%c", (headers[n] & (unsigned int) 0xFF000000) >> 24);
   }

   // These next 4 characters are for the biPlanes and biBitCount fields.

   fprintf(outfile, "%c", 1);
   fprintf(outfile, "%c", 0);
   fprintf(outfile, "%c", 24);
   fprintf(outfile, "%c", 0);

   for (n = 7; n <= 12; n++)
   {
      fprintf(outfile, "%c", headers[n] & 0x000000FF);
      fprintf(outfile, "%c", (headers[n] & 0x0000FF00) >> 8);
      fprintf(outfile, "%c", (headers[n] & 0x00FF0000) >> 16);
      fprintf(outfile, "%c", (headers[n] & (unsigned int) 0xFF000000) >> 24);
   }

   // Headers done, now write the data...

   for (y = height - 1; y >= 0; y--)  // BMPs are written bottom to top.
   {
      for (x = 0; x <= width - 1; x++)
      {
         // Also, it's written in (b,g,r) format...

         fprintf(outfile, "%c", pixels[(x * 3) + 2 + (y * width * 3)]);
         fprintf(outfile, "%c", pixels[(x * 3) + 1 + (y * width * 3)]);
         fprintf(outfile, "%c", pixels[(x * 3) + 0 + (y * width * 3)]);
      }
      if (extrabytes) // See above - BMP lines must be of lengths divisible by 4
      {
         for (n = 1; n <= extrabytes; n++)
         {
            fprintf(outfile, "%c", 0);
         }
      }
   }

   fclose (outfile);
   return;
}

Create a new file in your project and paste this into it.

The Epiphany code is:

#define set_red(n) (5*n  )
#define set_green(n) (20*n  )
#define set_blue(n) 0

__kernel void mandel_kern(
   int iterations,
   int width,
   float startx, 
   float starty, 
   float dx, 
   float dy, 
   __global uchar* pixels
)
{
int threeXwidth = 3 * width;
   unsigned char line[threeXwidth];
   int i, j, n, m, pixelBase;
   float x, y, r, s, nextr, nexts;

   j = get_global_id(0);

   for (i = 0; i < width; i++) 
   {

      x = startx + i*dx;
      y = starty + j*dy;

      r = x; 
      s = y;

      for (n = 0; n < iterations; n++) 
      {

         nextr = ((r * r) - (s * s)) + x;
         nexts = (2 * r * s) + y;
         r = nextr;
         s = nexts;
         
         if ((r * r) + (s * s) > 4) break;

      }

      if (n == iterations) n=0;

      line[(i * 3) + 0 ] = min(255,set_red(n));
      line[(i * 3) + 1 ] = min(255,set_green(n));
      line[(i * 3) + 2 ] = min(255,set_blue(n));   

   }

pixelBase = j * threeXwidth;
  for (m =0; m < threeXwidth; m++)
  pixels[pixelBase + m] = line[m];

}

Paste it into a new file called mandelbrot.cl. The compiler used by eclipse won't do anything with this file, it is here just for convenience.

1.4 Save and Build

If you've gotten this far you should now be able to save and build the program without any problems.

2 The Fun Starts Here


2.1 Debugging Host Code

Nothing works first time right? Right! I know that we are all geniuses but even we slip up occasionally.

Compile time on the host side is pretty much as usual. Get the includes and libraries right and it will all compile and link.

Run time debugging is when things start to get tricky. Prior to the execution of any code on the Epiphany accelerator using the Eclipse front end to the gdb debugger works as you would expect. However, the Epiphany expects the calling application to have root privileges. I started my project as an ordinary user (linaro) and when it came to starting Eclipse as root it got a bit sticky. First, logging in as root in a terminal window and then calling Eclipse didn't work at all and typing sudo eclipse from Run on the start menu got Eclipse running but then struck all sorts of permission problems. I reverted to old school debugging.

Notice I have opened a file called ./debug and I write a line to it pretty much before every call in the program. This is the cleaned up version of the debug writes that I used to figure out what was going on. I write them to a file to separate my output from the text that the stdcl libraries produce on the console (especially when the cl code actually runs). This is generally enough and if the Brown Deer documentation was a little more comprehensive would see you through.

Another lack in the documentation is how you go about checking if your call actually worked. (Notice that the only check I do is to open the cl file to see if it is there and readable.) If something does not work you get a Segmentation Fault a bit later. The Brown Deer's stdcl library takes away a lot of the verbosity of "generic" OpenCL code. For example, you don't have to create a context of your own, the library provides you with the stdacc global variable and that is the context for the Epiphany accelerator. However, what is lacking in the documentation is how to check if the calls worked e.g. if the clopen command found and successfully compiled the code. It might be there. If I find it I'll edit this post or write a new one.


2.2 Debugging Epiphany Code

I'm glad that I didn't really have to do this much.

Compile time with the JIT compiler is at run time for the host code. The output for the compiler is displayed onto the calling console. It reminded me of learning Pascal at university in the early 80's. I was fortunate that the cl code was correct enough for the compile time error to be obvious and with a little nutting out the changes I had to make were not too onerous. I think the lesson here is "Keep your kernels as simple as you can".

For run time debugging on the Epiphany there is a version of gdb called e-gdb. I did have not to go that far for this project but I think it will be key when things get more complicated. Andreas Olofsson got back to me about the developments in this area. There is significant work being done and things will get a lot easier in the near future. I think that for the moment, getting your code running as best you can in a friendly environment before you package it up into a kernel. Just make sure that you are only using the somewhat limited functions that are available if cl (the little copy loop at the end of the cl code above replaces a memcpy command in the original code).

I think I'll write another blog post once I get some idea about debugging on the Epiphany.


3. Final Thoughts


This guide is really only a toe in the water. So far I have only used the stdcl library and not gone anywhere near the COPRTHR library. There's lots of good stuff there that I have not gotten into yet.

The ARM cores on the Parallella are fine for a little snippet like the example but a bit too low powered for significant development. Cross compilation from a more powerful machine is the next step and I hope to write another blog entry on that in the near future.

Oh and one final thing... the output:



Friday, 20 December 2013

Sorting of Spatial Data for Meshed Multi-Core Processors

Sorting of Spatial Data for Meshed Multi-Core Processors

Many algorithms that operate on two dimensional data sets are well suited to multi-core processors such as the Parallella. The nature of raw two dimensional data is that it generally arrives as an unsorted “point cloud” where neither the x or y axis are sorted. In order to gain the most efficient use of the processor, the points must be sorted so that each core can work on a subset of the space while minimising the need for inter-core communication. Therefore it makes sense to start the process with an efficient, two dimensional sort where the data ends up evenly distributed amongst the cores in clusters of “near neighbours”.
The target of this and subsequent blogs will be a design for Delaunay Triangulation. The target architecture is the Parallella board. I believe that the Epiphany III and Epiphany IV chips on the Parallella is well suited to spatial processing given that its cores are organised in a mesh, i.e. each core is connected to it's nearest neighbours to the north, south, east and west via multiply high speed buses. Thus, if close data points are clustered on each core and neighbouring clusters are on neighbouring cores, the distance to be travelled when the clusters are finally matched up will be minimised.

Goals of the Algorithm

The goals of the sorting algorithm are:
  • The x value of each point in a cluster are less than those in the clusters on the cores to the east and greater than those on the cores to the west.
  • The y value of each point in a cluster are less than those in the clusters on the cores to the north and greater than those in the clusters on the cores to the south.
  • The points are evenly distributed between cores
In addition we should take care that no unnecessary steps are introduced so that we end up with the most efficient process over all.

Distributing the Points

Given that the point cloud is assumed to be random and we want to have the same number of points per core then a round robin distribution seems the best way to start. The can be done as the data is being read in which should give the cores enough time to complete at least some of the initial sort in parallel with the data load phase.


Initial Sort

Let's not reinvent the wheel here. I think that a quick sort of the x values into one index and the y values into another index is the way to go here.


Swap Sort

After the initial sort, each core will have an sorted group of points that could be from anywhere in the point cloud. The purpose of the swap sort is to push each point to the neighbouring cores where the point is with it's nearest neighbours, or at least closer to them. I'm using a push or write-only method of communication between cores because the write message is a lot quicker than the read message. The cores must swap points (i.e receive one point for every point sent) in order to preserve the even distribution of points.
The cores must start the process by passing their larges x and y values to the cores to the east and north respectively via the c-mesh (let's call this the MyLargestX message). If the lowest x value of the core to the east is smaller than the highest x value then the these two cores must initiate a swap. This can be done in response to the MyLargestX message (let's call the response MySmallesXY). Simultaneously, the cores can be swapping along the y axis with a MyLargestY message. The swap is then completed with a MyLargestXY call from the initial core passing both the x and y values to the points new “home”.
If the MyLargestX message contains an argument that is smaller than the smallest value on the receiving core then the receiving core does not initiate the swap. If, at a later time that core receives a point that is smaller than that received from the core to the south or west then the swap can be initiated.


End Conditions

The end conditions cannot be determined by any one core on behalf of all cores because each core only knows it's own values and the largest values of the cores to the south and west. Therefore, each core must indicate to the controlling application that it is has received a MyLargestX or MyLargestY that does not cause it to initiate a swap. This is only a tentative state because a core that is currently “happy” (i.e. has values that are larger than the cores to the south and west) may receive a point from the north or east that it must pass on. Therefore the controlling application must only terminate the process when all cores are “happy”.


Potential Problems

Sub-Optimal Swapping Pathways

Because each core can only make decisions based on the data that is currently in it's local memory there may be some cases that swaps are initiated early in the process are then undone due to some smaller values arriving from cores further to the north. Right now I can't see how this can be avoided.
Similarly, a swap to the west may be undone at a later stage after some swaps to the north. This could be avoided by swapping based on the x values (north-south) first and, when x is sorted, sorting on y (east-west). This would also free up some memory on each core given that there would only be one sort index needed at a time (or indeed the points could be moved around in memory removing the need for an index at all).


Skewed Data Sets

This kind of sorting will end up with the point cloud spread out into a square. This is fine if the original data set is also roughly square or at least regularly distributed into a regular sort of shape.




Diagram 1 – Near neighbours stay near.

This distribution probably will not work so well if the point cloud is not regular. Points that are not true near neighbours may end up on the same core and points that are near neighbours may end up on distant cores.



Diagram 2 – Distant points that are not near neighbours get squashed onto the same core(s)



Diagram 3 – Near neighbours get separated onto distant cores



In this case the user must be aware of the limitations of the sort and the algorithm that uses the sorted data must be robust enough to handle this situation.


Next: Using the sorted data

My next entry will a description of how to Triangulate the sorted points.
Don't hold your breath – the Parallella boards are being shipped and when I get mine, I'll be writing code, not english.