Computing SelfOrganizing Maps in a Massively Parallel Way with CUDA. Part 2: Algorithms
In the previous post I spoke briefly about motivations for implementing selforganizing maps in F# using GPU with CUDA. I have finally been able to outperform a single threaded C++ implementation by a factor of about 1.5. This is quite modest, but on the other hand rather impressive since we started out by being 60 times slower. At this point I am bidding farewell to F# and switching to C++. It will be interesting to see how this works out, but here are my initial algorithms.
So, we are parallelizing the following:
member this.GetBMU (node : Node) = let min = ref Double.MaxValue let minI = ref 1 let minJ = ref 1 this.somMap > Array2D.iteri (fun i j e > let dist = getDistance e node this.Metric if dist < !min then min := dist; minI := i; minJ := j) !minI, !minJ
Here somMap is just a twodimensional array of Nodes. And a Node is simply an array of float’s (double’s in C#), of the size equal to dimensionality of our space.
The code for getDistance is also simple:
let getDistanceEuclidian (x : Node) (y : Node) = Math.Sqrt([0..x.Dimension  1] > Seq.fold(fun sq i > sq + (x.[i]  y.[i]) ** 2.) 0.) let getDistance (x : Node) (y : Node) metric = if x.Dimension <> y.Dimension then failwith "Dimensions must match" else match metric with  Euclidian > getDistanceEuclidian x y  Taxicab > getDistanceTaxicab x y
In my parallel version I am only using Euclidian distance for simplicity.
Parallel Algorithms
As I have mentioned above, Alea.cuBase (my F#toCUDA framework) does not support multidimensional arrays (as well as shared memory or calling a kernel from a kernel), so this is the price I had to pay for sticking to my favorite language. Given these constraints, here is what I came up with:
Nodebynode
The very first idea (that proved also the very best) is quite intuitive. To find BMUs for the entire set of nodes, just take each individual node, find its BMU in parallel, repeat.
The map is first flattened into a single dimension, where if the cell coordinates were (i, j, k), they are mapped to (i * j * k + j * k + k). All distance subcalculations can be performed in one shot, and then one thread finishes them up. By “subcalculation” I mean computing (node(i) – map(j)) ** 2. These are then added up to the distance squared (I don’t calculate sqrt, since I don’t need it to find the minimum).
So, here is the implementation:
let pDistances = cuda { let! kernel = <@ fun nodeLen len (node : DevicePtr<float>) (map : DevicePtr<float>) (distances : DevicePtr<float>) (minDist : DevicePtr<float>) (minIndex : DevicePtr<int>) > // index into the original map, assuming // a node is a single entity let mapI = blockIdx.x * blockDim.x // actual position of the node component in the map let i = mapI + threadIdx.x if i < len then // index into the node let j = threadIdx.x % nodeLen distances.[i] < (map.[i]  node.[j]) * (map.[i]  node.[j]) if threadIdx.x = 0 then __syncthreads() let mutable thread = 0 minIndex.[blockIdx.x] < 1 // find the minimum among threads while mapI + thread < len && thread < blockDim.x do let k = mapI + thread let mutable sum = 0. for j = 0 to nodeLen  1 do sum < sum + distances.[k + j] if minDist.[blockIdx.x] > sum  minIndex.[blockIdx.x] < 0 then minDist.[blockIdx.x] < sum minIndex.[blockIdx.x] < k / nodeLen thread < thread + nodeLen @> > defineKernelFunc let diagnose (stats:KernelExecutionStats) = printfn "gpu timing: %10.3f ms %6.2f%% threads(%d) reg(%d) smem(%d)" stats.TimeSpan (stats.Occupancy * 100.0) stats.LaunchParam.BlockDim.Size stats.Kernel.NumRegs stats.Kernel.StaticSharedMemBytes return PFunc(fun (m:Module) (nodes : float [] list) (map : float []) > let kernel = kernel.Apply m let nodeLen = nodes.[0].Length let chunk = map.Length let nt = (256 / nodeLen) * nodeLen // number of threads divisible by nodeLen let nBlocks = (chunk + nt  1)/ nt //map.Length is a multiple of nodeLen by construction use dMap = m.Worker.Malloc(map) use dDist = m.Worker.Malloc<float>(map.Length) use dMinIndices = m.Worker.Malloc<int>(nBlocks * nodes.Length) use dMinDists = m.Worker.Malloc<float>(nBlocks * nodes.Length) use dNodes = m.Worker.Malloc(nodes.SelectMany(fun n > n :> float seq).ToArray()) let lp = LaunchParam(nBlocks, nt) //> Engine.setDiagnoser diagnose nodes > List.iteri (fun i node > kernel.Launch lp nodeLen chunk (dNodes.Ptr + i * nodeLen) dMap.Ptr dDist.Ptr (dMinDists.Ptr + i * nBlocks) (dMinIndices.Ptr + i * nBlocks)) let minDists = dMinDists.ToHost() let indices = dMinIndices.ToHost() let mins = (Array.zeroCreate nodes.Length) for i = 0 to nodes.Length  1 do let baseI = i * nBlocks let mutable min = minDists.[baseI] mins.[i] < indices.[baseI] for j = 1 to nBlocks  1 do if minDists.[baseI + j] < min then min <minDists.[baseI + j] mins.[i] < indices.[baseI + j] mins ) }
To get as much parallelism as possible, I start with 256 threads per block (maximum on Kepler is 1024, but 256 gives the best performance). Since each block of threads is going to compute as many distances as possible, I try to allocate the maximum number of threads divisible by node dimensionality. In my case: 256 / 12 * 256 = 252. Pretty good, we get almost an optimal number of threads per block.
The number of blocks are computed from the length of the map. I want to split this in an optimal way since each block is scheduled on one SP, and I have 192 of those I don’t want them to idle. The algorithm leans towards parallelizing calculations relative to the map (see picture above, I want all those “arrows” to be executed in parallel), so the number of blocks will be (somArray.length + nt – 1) / nt (nt is the number of threads – 252, somArray.length is 200 * 200 * 12 = 480000, the formula above takes into account the fact that this number may not be a multiple of 252, in which case we will need one more incomplete block). My block size is 1905 – pretty good, CUDA devs recommend that to be at least twice the number of multiprocessors. It is necessary to hide latency, which is killer in this development paradigm (you need to rely on your PCI slot to transfer data).
One weird thing here is that I have to allocate a relatively large dDist arrray. This is where temporary values of (map(i) – node(j))**2 go. In reality (if I were writing in C++) I would not do that. I would just use the superfast shared memory for this throwaway array. I could not get that to work with Alea, although the documentation says it is supported.
Another thing: the call to __syncthreads() is issued under a conditional statement. This, in general, is a horrible idea, because in the SIMT (single instruction multiple threads), threads march “in sync” so to speak, instruction by instruction. Thus, doing what I have done may lead to things hanging as some threads may take different branches and other threads will wait forever. Here, however, we are good, because the only way to go if the condition evaluates to false is to simply quit.
The kernel is called in a loop. One call for each node: 10000 passes (vs 10000 * 40000) in the sequential case. I also make sure that all of my memory is allocated once and everything I need is copied there. Not doing that leads to disastrous consequences, since all you would have in that case is pure latency.
The code is selfexplanatory. Once everyone has computed their part of the distance, these parts are assembled by the 0th thread of the block. That thread is responsible for calculating the final square of the distance, comparing the result of what we already have in the minDist array for this map node and storing that result.
Nodebynode optimized
And here lies the mistake that makes this algorithm lose out on performance: there is no need to delegate this work to the 0th thread (looping over all 252 threads of the block). It is enough to delegate that to each “threadIdx.x % nodeLen = 0″‘s thread of the block, so now 21 threads do this work in parallel, each for only 12 threads. Algorithm rewritten this way outperforms everything else I could come up with.
Here is the rewrite of the kernel function:
let! kernel = <@ fun nodeLen len (node : DevicePtr<float>) (map : DevicePtr<float>) (distances : DevicePtr<float>) (minDist : DevicePtr<float>) (minIndex : DevicePtr<int>) > // index into the original map, assuming // a node is a single entity let mapI = blockIdx.x * blockDim.x // actual position of the node component in the map let i = mapI + threadIdx.x if i < len then // index into the node let j = threadIdx.x % nodeLen distances.[i] < (map.[i]  node.[j]) * (map.[i]  node.[j]) __syncthreads() if threadIdx.x % nodeLen = 0 then minIndex.[blockIdx.x] < 1 // find the minimum among threads let k = mapI + threadIdx.x let mutable sum = 0. for j = k to k + nodeLen  1 do sum < sum + distances.[j] if minDist.[blockIdx.x] > sum  minIndex.[blockIdx.x] < 0 then minDist.[blockIdx.x] < sum minIndex.[blockIdx.x] < k / nodeLen @> > defineKernelFunc
This kernel function only stores “winning” (minimal) distances within each section of the map. Now they need to be reduced to one value per node. There are 1905 blocks and 10000 nodes. 10000 minimum values are computed in one pass over the 1905 * 10000length array of accumulated minimal distances:
let minDists = dMinDists.ToHost() let indices = dMinIndices.ToHost() let mins = (Array.zeroCreate nodes.Length) for i = 0 to nodes.Length  1 do let baseI = i * nBlocks let mutable min = minDists.[baseI] mins.[i] < indices.[baseI] for j = 1 to nBlocks  1 do if minDists.[baseI + j] < min then min <minDists.[baseI + j] mins.[i] < indices.[baseI + j] mins
And we are done. The improved version beats all the rest of my algorithms. Since all of these performance improvements looked so “juicy”, I decided it was high time to abandon managed code and go back to the C++ world. Especially since writing C++ code is not going to be so different: no annoying loops to write, just express it linearly and reap the massively parallel goodness.

October 6, 2013 at 1:02 pmF# Weekly #40, 2013  Sergey Tihon's Blog