<p>Here is some source code snip from our sample gallery you can find on the Alea GPU sample gallery <a href="http://www.quantalea.com/gallery" rel="nofollow">http://www.quantalea.com/gallery</a>.</p>

<p>Here is the single stage algorithm. It is not the fastest but reasonably simple to understand.</p>

<pre><code>public static class FloydWarshallSingleStage
{
    const int BlockWidth = 16;

    /// &lt;summary&gt;
    /// Kernel for parallel Floyd Warshall algorithm on GPU.
    /// &lt;/summary&gt;
    /// &lt;param name="u"&gt;Number vertex of which is performed relaxation paths [v1, v2]&lt;/param&gt;
    /// &lt;param name="n"&gt;Number of vertices in the graph G:=(V,E), n := |V(G)|&lt;/param&gt;
    /// &lt;param name="d"&gt;Matrix of shortest paths d(G)&lt;/param&gt;
    /// &lt;param name="p"&gt;Matrix of predecessors p(G)&lt;/param&gt;
    public static void KernelSingleStage(int u, int[,] d, int[,] p)
    {
        var n = d.GetLength(0);
        var v1 = blockDim.y * blockIdx.y + threadIdx.y;
        var v2 = blockDim.x * blockIdx.x + threadIdx.x;

        if (v1 &lt; n &amp;&amp; v2 &lt; n)
        {
            var newPath = d[v1, u] + d[u, v2];
            var oldPath = d[v1, v2];
            if (oldPath &gt; newPath)
            {
                d[v1, v2] = newPath;
                p[v1, v2] = p[u, v2];
            }
        }
    }

    [GpuManaged]
    public static void Run(Gpu gpu, int[,] d, int[,] p)
    {
        var n = d.GetLength(0);
        var gridDim = new dim3((n - 1) / BlockWidth + 1, (n - 1) / BlockWidth + 1, 1);
        var blockDim = new dim3(BlockWidth, BlockWidth, 1);
        var lp = new LaunchParam(gridDim, blockDim);
        for (var u = 0; u &lt; n; u++)
        {
            gpu.Launch(KernelSingleStage, lp, u, d, p);
        }
    }
}
</code></pre>

<p>The multi-stage is more complicated and the listing is longer. Here I paste a version that uses automatic memory management which simplifies the code quite a bit but also has some performance implication. The multi-stage version uses three kernels to complete the job and uses tiling to improve memory access.</p>

<pre><code>public class FloydWarshallMultiStage
{
    private const int None = -1;
    private const int Inf = 1061109567;

    //[GpuParam]
    //private readonly Constant&lt;int&gt; BlockSize;
    //[GpuParam]
    //private readonly Constant&lt;int&gt; ThreadSize;
    //[GpuParam]
    //private readonly Constant&lt;int&gt; VirtualBlockSize;

    private const int BlockSize = 16;
    private const int ThreadSize = 2;
    private const int VirtualBlockSize = BlockSize*ThreadSize;

    public FloydWarshallMultiStage(int blockSize, int threadSize)
    {
        //BlockSize = new Constant&lt;int&gt;(blockSize);
        //ThreadSize = new Constant&lt;int&gt;(threadSize);
        //VirtualBlockSize = new Constant&lt;int&gt;(blockSize * threadSize);
    }

    /// &lt;summary&gt;
    /// Kernel for parallel Floyd Warshall algorithm on GPU computing independent blocks.
    /// &lt;/summary&gt;
    /// &lt;param name="block"&gt;Number block of which is performed relaxation paths [v1, v2]&lt;/param&gt;
    /// &lt;param name="n"&gt;Number of vertices in the graph G:=(V,E), n := |V(G)|&lt;/param&gt;
    /// &lt;param name="pitch"&gt;Width to get to next row in number of int&lt;/param&gt;
    /// &lt;param name="d"&gt;Matrix of shortest paths d(G)&lt;/param&gt;
    /// &lt;param name="p"&gt;Matrix of predecessors p(G)&lt;/param&gt;
    public void KernelPhaseOne(int block, int n, int pitch, int[,] d, int[,] p)
    {
        var newPred = 0;

        var tx = threadIdx.x;
        var ty = threadIdx.y;

        var v1 = VirtualBlockSize*block + ty;
        var v2 = VirtualBlockSize*block + tx;

        var primaryD = __shared__.Array2D&lt;int&gt;(VirtualBlockSize, VirtualBlockSize);
        var primaryP = __shared__.Array2D&lt;int&gt;(VirtualBlockSize, VirtualBlockSize);

        if (v1 &lt; n &amp;&amp; v2 &lt; n)
        {
            primaryD[ty, tx] = d[v1, v2];
            primaryP[ty, tx] = p[v1, v2];
            newPred = primaryP[ty, tx];
        }
        else
        {
            primaryD[ty, tx] = Inf;
            primaryP[ty, tx] = None;
        }

        DeviceFunction.SyncThreads();

        for (var i = 0; i &lt; VirtualBlockSize; i++)
        {
            var newPath = primaryD[ty, i] + primaryD[i, tx];

            DeviceFunction.SyncThreads();

            if (newPath &lt; primaryD[ty, tx])
            {
                primaryD[ty, tx] = newPath;
                newPred = primaryP[i, tx];
            }

            DeviceFunction.SyncThreads();

            primaryP[ty, tx] = newPred;
        }

        if (v1 &lt; n &amp;&amp; v2 &lt; n)
        {
            d[v1, v2] = primaryD[ty, tx];
            p[v1, v2] = primaryP[ty, tx];
        }
    }

    /// &lt;summary&gt;
    /// Kernel for parallel Floyd Warshall algorithm on GPU to compute block depending on a single independent block.
    /// &lt;/summary&gt;
    /// &lt;param name="block"&gt;Number block of which is performed relaxation paths [v1, v2]&lt;/param&gt;
    /// &lt;param name="n"&gt;Number of vertices in the graph G:=(V,E), n := |V(G)|&lt;/param&gt;
    /// &lt;param name="pitch"&gt;&lt;/param&gt;
    /// &lt;param name="d"&gt;Matrix of shortest paths d(G)&lt;/param&gt;
    /// &lt;param name="p"&gt;Matrix of predecessors p(G)&lt;/param&gt;
    public void KernelPhaseTwo(int block, int n, int pitch, int[,] d, int[,] p)
    {
        if (blockIdx.x == block) return;
        var newPath = 0;
        var newPred = 0;

        var tx = threadIdx.x;
        var ty = threadIdx.y;

        var v1 = VirtualBlockSize*block + ty;
        var v2 = VirtualBlockSize*block + tx;

        var primaryD = __shared__.Array2D&lt;int&gt;(VirtualBlockSize, VirtualBlockSize);
        var currentD = __shared__.Array2D&lt;int&gt;(VirtualBlockSize, VirtualBlockSize);
        var primaryP = __shared__.Array2D&lt;int&gt;(VirtualBlockSize, VirtualBlockSize);
        var currentP = __shared__.Array2D&lt;int&gt;(VirtualBlockSize, VirtualBlockSize);

        if (v1 &lt; n &amp;&amp; v2 &lt; n)
        {
            primaryD[ty, tx] = d[v1, v2];
            primaryP[ty, tx] = p[v1, v2];
        }
        else
        {
            primaryD[ty, tx] = Inf;
            primaryP[ty, tx] = None;
        }

        // load i-aligned singly dependent blocks
        if (blockIdx.y == 0)
        {
            v1 = VirtualBlockSize*block + ty;
            v2 = VirtualBlockSize*blockIdx.x + tx;
        }
        // load j-aligned singly dependent blocks
        else
        {
            v1 = VirtualBlockSize*blockIdx.x + ty;
            v2 = VirtualBlockSize*block + tx;
        }

        if (v1 &lt; n &amp;&amp; v2 &lt; n)
        {
            currentD[ty, tx] = d[v1, v2];
            currentP[ty, tx] = p[v1, v2];
            newPred = currentP[ty, tx];
        }
        else
        {
            currentD[ty, tx] = Inf;
            currentP[ty, tx] = None;
        }

        DeviceFunction.SyncThreads();

        // compute i-aligned singly dependent blocks
        if (blockIdx.y == 0)
        {
            for (var i = 0; i &lt; VirtualBlockSize; i++)
            {
                newPath = primaryD[ty, i] + currentD[i, tx];

                DeviceFunction.SyncThreads();

                if (newPath &lt; currentD[ty, tx])
                {
                    currentD[ty, tx] = newPath;
                    newPred = currentP[i, tx];
                }

                DeviceFunction.SyncThreads();

                currentP[ty, tx] = newPred;
            }
        }
        // compute j-aligned singly dependent blocks
        else
        {
            for (var i = 0; i &lt; VirtualBlockSize; i++)
            {
                newPath = currentD[ty, i] + primaryD[i, tx];

                DeviceFunction.SyncThreads();
                if (newPath &lt; currentD[ty, tx])
                {
                    currentD[ty, tx] = newPath;
                    currentP[ty, tx] = primaryP[i, tx];
                }

                DeviceFunction.SyncThreads();
            }
        }

        if (v1 &lt; n &amp;&amp; v2 &lt; n)
        {
            d[v1, v2] = currentD[ty, tx];
            p[v1, v2] = currentP[ty, tx];
        }
    }

    /// &lt;summary&gt;
    /// Kernel for parallel Floyd Warshall algorithm on GPU to compute dependent block depending on the singly dependent blocks.
    /// &lt;/summary&gt;
    /// &lt;param name="block"&gt;Number block of which is performed relaxation paths [v1, v2]&lt;/param&gt;
    /// &lt;param name="n"&gt;Number of vertices in the graph G:=(V,E), n := |V(G)|&lt;/param&gt;
    /// &lt;param name="pitch"&gt;&lt;/param&gt;
    /// &lt;param name="d"&gt;Matrix of shortest paths d(G)&lt;/param&gt;
    /// &lt;param name="p"&gt;Matrix of predecessors p(G)&lt;/param&gt;
    public void KernelPhaseThree(int block, int n, int pitch, int[,] d, int[,] p)
    {
        if (blockIdx.x == block || blockIdx.y == block) return;

        var tx = threadIdx.x*ThreadSize;
        var ty = threadIdx.y*ThreadSize;

        var v1 = blockDim.y*blockIdx.y*ThreadSize + ty;
        var v2 = blockDim.x*blockIdx.x*ThreadSize + tx;

        var primaryRowD = __shared__.Array2D&lt;int&gt;(BlockSize*ThreadSize, BlockSize*ThreadSize);
        var primaryColD = __shared__.Array2D&lt;int&gt;(BlockSize*ThreadSize, BlockSize*ThreadSize);
        var primaryRowP = __shared__.Array2D&lt;int&gt;(BlockSize*ThreadSize, BlockSize*ThreadSize);

        var v1Row = BlockSize*block*ThreadSize + ty;
        var v2Col = BlockSize*block*ThreadSize + tx;

        // load data for virtual block
        for (var i = 0; i &lt; ThreadSize; i++)
        {
            for (var j = 0; j &lt; ThreadSize; j++)
            {
                var idx = tx + j;
                var idy = ty + i;

                if (v1Row + i &lt; n &amp;&amp; v2 + j &lt; n)
                {
                    primaryRowD[idy, idx] = d[v1Row + i, v2 + j];
                    primaryRowP[idy, idx] = p[v1Row + i, v2 + j];
                }
                else
                {
                    primaryRowD[idy, idx] = Inf;
                    primaryRowP[idy, idx] = None;
                }

                if (v1 + i &lt; n &amp;&amp; v2Col + j &lt; n)
                {
                    primaryColD[idy, idx] = d[v1 + i, v2Col + j];
                }
                else
                {
                    primaryColD[idy, idx] = Inf;
                }
            }
        }

        DeviceFunction.SyncThreads();

        // compute data for virtual block
        for (var i = 0; i &lt; ThreadSize; i++)
        {
            for (var j = 0; j &lt; ThreadSize; j++)
            {
                if (v1 + i &lt; n &amp;&amp; v2 + j &lt; n)
                {
                    var path = d[v1 + i, v2 + j];
                    var predecessor = p[v1 + i, v2 + j];

                    var idy = ty + i;
                    var idx = tx + j;

                    for (var k = 0; k &lt; BlockSize*ThreadSize; k++)
                    {
                        var newPath = primaryColD[idy, k] + primaryRowD[k, idx];
                        if (path &gt; newPath)
                        {
                            path = newPath;
                            predecessor = primaryRowP[k, idx];
                        }
                    }
                    d[v1 + i, v2 + j] = path;
                    p[v1 + i, v2 + j] = predecessor;
                }
            }
        }
    }

    /// &lt;summary&gt;
    /// Parallel multi-stage Floyd Warshall algorithm on GPU.
    /// &lt;/summary&gt;
    /// &lt;param name="gpu"&gt;The GPU on which the kernels should run&lt;/param&gt;
    /// &lt;param name="n"&gt;Number of vertices in the graph G:=(V,E), n := |V(G)|&lt;/param&gt;
    /// &lt;param name="g"&gt;The graph G:=(V,E)&lt;/param&gt;
    /// &lt;param name="d"&gt;Matrix of shortest paths d(G)&lt;/param&gt;
    /// &lt;param name="p"&gt;Matrix of predecessors p(G)&lt;/param&gt;
    public void Run(Gpu gpu, int[,] d, int[,] p, bool verbose = false)
    {
        var n = d.GetLength(0);
        var gridDim1 = new dim3(1, 1, 1);
        var gridDim2 = new dim3((n - 1)/VirtualBlockSize + 1, 2, 1);
        var gridDim3 = new dim3((n - 1)/VirtualBlockSize + 1, (n - 1)/VirtualBlockSize + 1, 1);

        var blockDim1 = new dim3(VirtualBlockSize, VirtualBlockSize, 1);
        var blockDim2 = new dim3(VirtualBlockSize, VirtualBlockSize, 1);
        var blockDim3 = new dim3(BlockSize, BlockSize, 1);

        var numOfBlock = (n - 1)/VirtualBlockSize + 1;
        var pitchInt = n;

        if (verbose)
        {
            Console.WriteLine($"|V| {n}");
            Console.WriteLine($"Phase 1: grid dim {gridDim1} block dim {blockDim1}");
            Console.WriteLine($"Phase 2: grid dim {gridDim2} block dim {blockDim2}");
            Console.WriteLine($"Phase 3: grid dim {gridDim3} block dim {blockDim3}");
        }

        for (var block = 0; block &lt; numOfBlock; block++)
        {
            gpu.Launch(KernelPhaseOne, new LaunchParam(gridDim1, blockDim1), block, n, pitchInt, d, p);
            gpu.Launch(KernelPhaseTwo, new LaunchParam(gridDim2, blockDim2), block, n, pitchInt, d, p);
            gpu.Launch(KernelPhaseThree, new LaunchParam(gridDim3, blockDim3), block, n, pitchInt, d, p);
        }
    }
}
</code></pre>

<p>The version with explicit memory management you better download the sample from <a href="http://www.quantalea.com/gallery" rel="nofollow">http://www.quantalea.com/gallery</a> and search for Floyd-Warshall. </p>

<p>Hope that answers the question. </p>

<p>The implementation is based on the following paper:</p>

<p>Ben Lund, Justin W Smith, A Multi-Stage CUDA Kernel for Floyd-Warshall, 2010.</p>

<p><a href="https://arxiv.org/abs/1001.4108" rel="nofollow">https://arxiv.org/abs/1001.4108</a></p>
Floyd Warshall using Alea GPU
<p>I've been trying to use Alea GPU to write the parallel Floyd-Warshall algorithm in F#, and basing myself on the CUDA code another user presented here</p>

<p><a href="https://stackoverflow.com/questions/19861532/the-floyd-warshall-algorithm-in-cuda">The Floyd-Warshall algorithm in CUDA</a></p>

<p>I wrote the following simple implementation</p>

<pre><code>type FWModule&lt;'T&gt;(target:GPUModuleTarget, tileDim:int) =
inherit GPUModule(target)

    [&lt;Kernel;ReflectedDefinition&gt;]
    member this.FloydWKernel (width:int) (k:int) (data:deviceptr&lt;float&gt;) =

        let col = blockIdx.x * blockDim.x + threadIdx.x
        let row = blockIdx.y

        if col &gt;= width then ()  //out of bounds
        let index  = width * row + col
        let best = __shared__.Variable&lt;float&gt;()

        if threadIdx.x = 0 then best := data.[width*row+k]
        __syncthreads()

        let tmp = data.[k*width+col]

        let candidate = !best + tmp
        data.[index] &lt;- min data.[index] candidate

    member this.LaunchParams width =
        let blockdim = dim3(tileDim)
        let griddim = dim3(divup width tileDim, width)
        LaunchParam(griddim, blockdim)

    member this.FloydW (width:int) (k:int) (data:deviceptr&lt;float&gt;) =
        let lp = this.LaunchParams width
        this.GPULaunch &lt;@ this.FloydWKernel @&gt; lp width k idata odata

    member this.FW(size:int, A:float[])=
        use deviceArr = this.GPUWorker.Malloc(A)
        for k in 0 .. size-1 do
           this.FloydW size k deviceArr.Ptr deviceArr.Ptr
        deviceArr.Gather()

let tileDim = 256

let apsp = new FWModule&lt;float&gt;(GPUModuleTarget.DefaultWorker, tileDim)
</code></pre>

<p>However, when the following lines are ran in <code>fsi</code></p>

<pre><code>let m = [|0.0     ; 5.0     ; 9.0     ; infinity;
          infinity; 0.0     ; 1.0     ; infinity;
          infinity; infinity; 0.0     ; 2.0;
          infinity; 3.0     ; infinity; 0.0|];;

apsp.FW (4,m);;
</code></pre>

<p>The output is</p>

<pre><code>[|0.0; 5.0; 6.0; 8.0;
  4.0; 0.0; 1.0; 3.0;
  3.0; 3.0; 0.0; 1.0;
  1.0; 1.0; 1.0; 0.0|]
</code></pre>

<p>Which it should not be given that the usual iterative, sequential <code>floydwarshall</code></p>

<pre><code>let floydwarshall (l:int, mat:float[]) =
   let a = Array.copy mat
   for k in 0 .. (l-1) do
      for i in 0 .. (l-1) do
         for j in 0 .. (l-1) do
            a.[i*l+j] &lt;- min a.[i*l+j] (a.[i*l+k] + a.[k*l+j])
   a
</code></pre>

<p>gives me</p>

<pre><code>floydwarshall (4,m);;
[|0.0     ; 5.0; 6.0; 8.0;
  infinity; 0.0; 1.0; 3.0;
  infinity; 5.0; 0.0; 2.0;
  infinity; 3.0; 4.0; 0.0|]
</code></pre>

<p>My question is, what's happening? </p>