A simple experiment to understand CUDAFunctionLoad

Question

In an attempt to answer my own question (Beginner CUDA: two dimensional blocks and two dimensional threads?).

Please excuse the long post, it contains for examples I constructed as an experiment to understand what is going on. (postscript: on the advice of @rojo, I am clarifying my question; so the post is getting longer. New edits in italics).

I've constructed the simplest CUDA code example I can think of, but I don't understand the results. I would be grateful if someone with more CUDAFunctionLoad experience than I could give me an explanation.

src = "__global__ void query(mint * idx, mint * tinx, mint * tiny, \
mint * binx, mint *biny, mint * bdimx, mint *bdimy, mint * gdimx, \
mint * gdimy, mint length)
   {
      int index = threadIdx.x + blockIdx.x*blockDim.x;
      if (index < length)
          {
          idx[index] = index;
          tinx[index] = threadIdx.x  ;
          tiny[index] = threadIdx.y;
          binx[index] = blockIdx.x  ;
          biny[index] = blockIdx.y;
          bdimx[index] = blockDim.x; 
          bdimy[index] = blockDim.y;
          gdimx[index] = gridDim.x; 
          gdimy[index] = gridDim.y;
          }
  }";

Use CUDAFunctionLoad with a {16,8} blockdim argument:

simpleCUDA = 
 CUDAFunctionLoad[src, 
  "query", {{_Integer, "Output"}, {_Integer, "Output"}, {_Integer, 
    "Output"}, {_Integer, "Output"}, {_Integer, "Output"}, {_Integer, 
    "Output"}, {_Integer, "Output"}, {_Integer, "Output"}, {_Integer, 
    "Output"}, _Integer}, {16, 8}, "ShellOutputFunction" -> Print]


headings = {"index", "thread index x", "thread index y", 
   "block index x", "block index y", "block dimension x", 
   "block dimension y", "grid dimension x", "grid dimension y"};
avec = Range[20];

Call this function without the optional threads last argument:

(Why are the grid dimensions {2,1}? The thread index y changes at index 16 as expected, but why are the indices 7 and then 3--are these chosen at random?)

MatrixForm@
 Transpose@
  Prepend[Transpose@
    simpleCUDA[avec, avec, avec, avec, avec, avec, avec, avec, avec, 
     20], headings]

With result:

Call the function with the optional threads final argument (e.g. 6):

(Why should the grid dimensions change at index 16? How does one understand the effect of the final argument 6?)

MatrixForm@
 Transpose@
  Prepend[Transpose@
    simpleCUDA[avec, avec, avec, avec, avec, avec, avec, avec, avec, 
     20, 6], headings]

Result:

Use CUDAFunctionLoad with a single blockdim argument (e.g., 8)

simpleCUDAalt = 
  CUDAFunctionLoad[{getSource["simple.cu"]}, 
   "query", {{_Integer, "Output"}, {_Integer, "Output"}, {_Integer, 
     "Output"}, {_Integer, "Output"}, {_Integer, "Output"}, {_Integer,
      "Output"}, {_Integer, "Output"}, {_Integer, 
     "Output"}, {_Integer, "Output"}, _Integer}, 8, 
   "ShellOutputFunction" -> Print];

Call this function without the optional thread argument:

MatrixForm@
 Transpose@
  Prepend[Transpose@
    simpleCUDAalt[avec, avec, avec, avec, avec, avec, avec, avec, 
     avec, 20], headings]

Result:

(Why are the grid dimensions now {1,1} when previously they were {2,1}? Why do the grid dimensions change at index 16? Is it that I am accessing memory that hasn't been cleared?)

It looks like I am only allowed to put 2 images in a single post--the editor is not allowing me to put in a 3rd. So, I will cut and paste Output

{
 {"index", 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 17, 
  18, 19, 20},
 {"thread index x", 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 
  15, 17, 18, 19, 20},
 {"thread index y", 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 
  17, 18, 19, 20},
 {"block index x", 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 17,
   18, 19, 20},
 {"block index y", 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 17,
   18, 19, 20},
 {"block dimension x", 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16,
   16, 16, 16, 16, 17, 18, 19, 20},
 {"block dimension y", 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
   17, 18, 19, 20},
 {"grid dimension x", 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
  17, 18, 19, 20},
 {"grid dimension y", 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
  17, 18, 19, 20}
}

Lastly, call the function with the additional threads argument (e.g., 12).

MatrixForm@
 Transpose@
  Prepend[Transpose@
    simpleCUDAalt[avec, avec, avec, avec, avec, avec, avec, avec, 
     avec, 20, 12], headings]

Result:

(Why do the block dimensions change as well here?)

{
 {"index", 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 17, 
  18, 19, 20},
 {"thread index x", 0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7, 
  17, 18, 19, 20},
 {"thread index y", 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
  17, 18, 19, 20},
 {"block index x", 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 17,
   18, 19, 20},
 {"block index y", 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 17,
   18, 19, 20},
 {"block dimension x", 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
   17, 18, 19, 20},
 {"block dimension y", 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
   17, 18, 19, 20},
 {"grid dimension x", 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 
  17, 18, 19, 20},
 {"grid dimension y", 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
  17, 18, 19, 20}
}

Answer 1

I think I've figured this block and grid dimension stuff out. First, here's a simpler kernel for checking:

dementiaSource = "
  __global__ void dementia( mint *out )
  {
    out[0] = blockDim.x;
    out[1] = blockDim.y;
    out[2] = blockDim.z;
    out[3] = gridDim.x;
    out[4] = gridDim.y;
    out[5] = gridDim.z;
  }";

This, of course, is highly redundant, with every thread poking the same values into the output array, but it does the job.

Make an instance with diagnostic block dimensions:

dementia = CUDAFunctionLoad[dementiaSource, "dementia", 
  {{_Integer, _, "Output"}}, {3, 5, 7}]

Now, invoke it with diagnostic dimensions:

dementia[ {0, 0, 0, 0, 0, 0}, {11, 13}]
(* {{3, 5, 7, 4, 3, 1}} *)

Note the the grid X dimension, 4, is Ceiling[11/3] and the grid Y dimension, 3, is Ceiling[13/5]. More experiments confirm that this appears to be a general rule.

My interpretation of this is that while CUDA's model of threading is six dimensional, CUDALink's model is basically two dimensional. CUDALink expects you to tile a two dimensional space with two dimensional blocks. The dimensions you provide to the CUDAFunction are not the grid dimensions, but the dimensions of the space you are tiling in elements (pixels). CUDALink rounds up the X and Y grid dimensions as needed to completely cover the space. The Z block dimension given CUDAFunctionLoad is used, but there is no corresponding Z dimension in the CUDAFunction parameters. To recover a five dimensional model, the dimensions provided to the CUDAFunction should be {gx*bx,gy*by} where gx is the X grid dimension, and bx is the X block dimension, etc. There appears to be no way to specify a grid Z dimension.