openCL reduction, and passing 2d array

Question

Here is the loop I want to convert to openCL.

for(n=0; n < LargeNumber; ++n) {    
    for (n2=0; n2< SmallNumber; ++n2) {
        A[n]+=B[n2][n];
    }                                                         
    Re+=A[n];       
}

And here is what I have so far, although, I know it is not correct and missing some things.

__kernel void openCL_Kernel( __global  int *A,
                         __global  int **B,  
                         __global  int *C, 
                         __global _int64 Re,
                                   int D) 
{

int i=get_global_id(0);
int ii=get_global_id(1);

A[i]+=B[ii][i];

//barrier(..); ?

Re+=A[i];

}

I'm a complete beginner to this type of thing. First of all I know that I can't pass a global double pointer to an openCL kernel. If you can, wait a few days or so before posting the solution, I want to figure this out for myself, but if you can help point me in the right direction I would be grateful.

Answer 1

Concerning your problem with passing doublepointers: That kind of problem is typically solved by copying the whole matrix (or whatever you are working on) into one continous block of memory and, if the blocks have different lengths passing another array, which contains the offsets for the individual rows ( so your access would look something like B[index[ii]+i]).

Now for your reduction down to Re: since you didn't mention what kind of device you are working on I'm going to assume its GPU. In that case I would avoid doing the reduction in the same kernel, since its going to be slow as hell the way you posted it (you would have to serialize the access to Re over thousands of threads (and the access to A[i] too). Instead I would write want kernel, which sums all B[*][i] into A[i] and put the reduction from A into Re in another kernel and do it in several steps, that is you use a reduction kernel which operates on n element and reduces them to something like n / 16 (or any other number). Then you iteratively call that kernel until you are down to one element, which is your result (I'm making this description intentionally vague, since you said you wanted to figure thinks out yourself).

As a sidenote: You realize that the original code doesn't exactly have a nice memory access pattern? Assuming B is relatively large (and much larger then A due to the second dimension) having the inner loop iterate over the outer index is going to create a lot of cachemisses. This is even worse when porting to the gpu, which is very sensitive about coherent memory access

So reordering it like this may massively increase performance:

for (n2=0; n2< SmallNumber; ++n2)
  for(n=0; n < LargeNumber; ++n)    
    A[n]+=B[n2][n];
for(n=0; n < LargeNumber; ++n)                                                 
  Re+=A[n];       

This is particulary true if you have a compiler which is good at autovectorization, since it might be able to vectorize that construct, but it's very unlikely to be able to do so for the original code (and if it can't prove that A and B[n2] can't refer to the same memory it can't turn the original code into this).