add: msbfs levels

Author: artemiipatov

src/GraphBLAS-sharp.Backend/Algorithms/MSBFS.fs

@@ -0,0 +1,136 @@
+namespace GraphBLAS.FSharp.Backend.Algorithms
+
+open Brahma.FSharp
+open FSharp.Quotations
+open GraphBLAS.FSharp
+open GraphBLAS.FSharp.Objects
+open GraphBLAS.FSharp.Objects.ClMatrix
+open GraphBLAS.FSharp.Objects.ArraysExtensions
+open GraphBLAS.FSharp.Objects.ClContextExtensions
+open GraphBLAS.FSharp.Backend.Matrix.LIL
+open GraphBLAS.FSharp.Backend.Matrix.COO
+
+module internal MSBFS =
+    module Levels =
+        let private updateFront (clContext: ClContext) workGroupSize =
+
+            let excludeValues = ClArray.excludeElements clContext workGroupSize
+
+            let excludeIndices = ClArray.excludeElements clContext workGroupSize
+
+            fun (queue: MailboxProcessor<_>) allocationMode (front: ClMatrix.COO<_>) (intersection: ClArray<int>) ->
+
+                let newRows = excludeIndices queue allocationMode intersection front.Rows
+
+                let newColumns = excludeIndices queue allocationMode intersection front.Columns
+
+                let newValues = excludeValues queue allocationMode intersection front.Values
+
+                match newRows, newColumns, newValues with
+                | Some rows, Some columns, Some values ->
+                    { Context = clContext
+                      Rows = rows
+                      Columns = columns
+                      Values = values
+                      RowCount = front.RowCount
+                      ColumnCount = front.ColumnCount }
+                    |> Some
+                | _ -> None
+
+        let private updateFrontAndLevels (clContext: ClContext) workGroupSize =
+
+            let updateFront = updateFront clContext workGroupSize
+
+            let mergeDisjoint = Matrix.mergeDisjoint clContext workGroupSize
+
+            let findIntersection = Intersect.findKeysIntersection clContext workGroupSize
+
+            fun (queue: MailboxProcessor<_>) allocationMode (levels: ClMatrix.COO<_>) (front: ClMatrix.COO<_>) ->
+
+                // Find intersection of levels and front indices.
+                let intersection = findIntersection queue DeviceOnly front levels
+
+                // Remove mutual elements
+                let newFront = updateFront queue allocationMode front intersection
+
+                intersection.Free queue
+
+                match newFront with
+                | Some f ->
+                    // Update levels
+                    let newLevels = mergeDisjoint queue levels f
+                    newLevels, newFront
+                | _ -> levels, None
+
+        let run<'a when 'a: struct>
+            (add: Expr<int -> int -> int option>)
+            (mul: Expr<'a -> int -> int option>)
+            (clContext: ClContext)
+            workGroupSize
+            =
+
+            let spGeMM =
+                Operations.SpGeMM.expand add mul clContext workGroupSize
+
+            let toCSRInPlace = Matrix.toCSRInPlace clContext workGroupSize
+
+            let updateFrontAndLevels = updateFrontAndLevels clContext workGroupSize
+
+            fun (queue: MailboxProcessor<Msg>) (matrix: ClMatrix<'a>) (source: int list) ->
+                let vertexCount = matrix.RowCount
+                let sourceVertexCount = source.Length
+
+                let startMatrix =
+                    source
+                    |> List.mapi (fun i vertex -> i, vertex, 1)
+
+                let mutable levels =
+                    startMatrix
+                    |> Matrix.ofList clContext DeviceOnly sourceVertexCount vertexCount
+
+                let mutable front =
+                    startMatrix
+                    |> Matrix.ofList clContext DeviceOnly sourceVertexCount vertexCount
+                    |> toCSRInPlace queue DeviceOnly
+
+                let mutable level = 0
+                let mutable stop = false
+
+                while not stop do
+                    level <- level + 1
+
+                    //Getting new frontier
+                    match spGeMM queue DeviceOnly matrix (ClMatrix.CSR front) with
+                    | None ->
+                        front.Dispose queue
+                        stop <- true
+                    | Some newFrontier ->
+                        front.Dispose queue
+                        //Filtering visited vertices
+                        match updateFrontAndLevels queue DeviceOnly levels newFrontier with
+                        | l, Some f ->
+                            front <- toCSRInPlace queue DeviceOnly f
+                            levels.Dispose queue
+                            levels <- l
+                            newFrontier.Dispose queue
+                        | _, None ->
+                            stop <- true
+                            newFrontier.Dispose queue
+
+                levels
+
+        let runSingleSourceMultipleTimes<'a when 'a: struct>
+            (add: Expr<int option -> int option -> int option>)
+            (mul: Expr<'a option -> int option -> int option>)
+            (clContext: ClContext)
+            workGroupSize
+            =
+
+            let SSBFS = BFS.singleSourceSparse add mul clContext workGroupSize
+
+            fun (queue: MailboxProcessor<Msg>) (matrix: ClMatrix<'a>) (source: int list) ->
+                source
+                |> List.map (SSBFS queue matrix)
+
+    module Parents =
+        let run = 0

src/GraphBLAS-sharp.Backend/Common/ClArray.fs

@@ -121,6 +121,37 @@ module ClArray =
 
             outputArray
 
+    /// <summary>
+    /// Copies all elements from source to destination array.
+    /// </summary>
+    /// <param name="clContext">OpenCL context.</param>
+    /// <param name="workGroupSize">Should be a power of 2 and greater than 1.</param>
+    let copyTo (clContext: ClContext) workGroupSize =
+        let copy =
+            <@ fun (ndRange: Range1D) (source: ClArray<'a>) (destination: ClArray<'a>) inputArrayLength ->
+
+                let i = ndRange.GlobalID0
+
+                if i < inputArrayLength then
+                    source.[i] <- destination.[i] @>
+
+        let program = clContext.Compile(copy)
+
+        fun (processor: MailboxProcessor<_>) (source: ClArray<'a>) (destination: ClArray<'a>) ->
+            if source.Length <> destination.Length then
+                failwith "The source array length differs from the destination array length."
+
+            let ndRange =
+                Range1D.CreateValid(source.Length, workGroupSize)
+
+            let kernel = program.GetKernel()
+
+            processor.Post(
+                Msg.MsgSetArguments(fun () -> kernel.KernelFunc ndRange source destination source.Length)
+            )
+
+            processor.Post(Msg.CreateRunMsg<_, _> kernel)
+
     /// <summary>
     /// Creates an array of the given size by replicating the values of the given initial array.
     /// </summary>
@@ -781,3 +812,83 @@ module ClArray =
             bitmap.Free processor
 
             result
+
+    /// <summary>
+    /// Builds a new array whose elements are the results of applying the given function
+    /// to each of the elements of the array.
+    /// </summary>
+    /// <param name="op">The function to transform elements of the array.</param>
+    /// <param name="clContext">OpenCL context.</param>
+    /// <param name="workGroupSize">Should be a power of 2 and greater than 1.</param>
+    let map<'a, 'b> (op: Expr<'a -> 'b>) (clContext: ClContext) workGroupSize = Map.map op clContext workGroupSize
+
+    /// <summary>
+    /// Builds a new array whose elements are the results of applying the given function
+    /// to each of the elements of the array.
+    /// </summary>
+    /// <param name="op">The function to transform elements of the array.</param>
+    /// <param name="clContext">OpenCL context.</param>
+    /// <param name="workGroupSize">Should be a power of 2 and greater than 1.</param>
+    let mapInPlace<'a> (op: Expr<'a -> 'a>) (clContext: ClContext) workGroupSize = Map.mapInPlace op clContext workGroupSize
+
+    /// <summary>
+    /// Builds a new array whose elements are the results of applying the given function
+    /// to the corresponding pairs of values, where the first element of pair is from the given array
+    /// and the second element is the given value.
+    /// </summary>
+    /// <param name="op">The function to transform elements of the array.</param>
+    /// <param name="clContext">OpenCL context.</param>
+    /// <param name="workGroupSize">Should be a power of 2 and greater than 1.</param>
+    let mapWithValue<'a, 'b, 'c> (clContext: ClContext) workGroupSize (op: Expr<'a -> 'b -> 'c>) = Map.mapWithValue clContext workGroupSize op
+
+    /// <summary>
+    /// Builds a new array whose elements are the results of applying the given function
+    /// to the corresponding elements of the two given arrays pairwise.
+    /// </summary>
+    /// <remarks>
+    /// The two input arrays must have the same lengths.
+    /// </remarks>
+    /// <param name="map">The function to transform the pairs of the input elements.</param>
+    /// <param name="clContext">OpenCL context.</param>
+    /// <param name="workGroupSize">Should be a power of 2 and greater than 1.</param>
+    let map2<'a, 'b, 'c> (map: Expr<'a -> 'b -> 'c>) (clContext: ClContext) workGroupSize = Map.map2 map clContext workGroupSize
+
+    /// <summary>
+    /// Fills the third given array with the results of applying the given function
+    /// to the corresponding elements of the first two given arrays pairwise.
+    /// </summary>
+    /// <remarks>
+    /// The first two input arrays must have the same lengths.
+    /// </remarks>
+    /// <param name="map">The function to transform the pairs of the input elements.</param>
+    /// <param name="clContext">OpenCL context.</param>
+    /// <param name="workGroupSize">Should be a power of 2 and greater than 1.</param>
+    let map2InPlace<'a, 'b, 'c> (map: Expr<'a -> 'b -> 'c>) (clContext: ClContext) workGroupSize = Map.map2InPlace map clContext workGroupSize
+
+    /// <summary>
+    /// Excludes elements, pointed by the bitmap.
+    /// </summary>
+    /// <param name="clContext">OpenCL context.</param>
+    /// <param name="workGroupSize">Should be a power of 2 and greater than 1.</param>
+    let excludeElements (clContext: ClContext) workGroupSize =
+
+        let invert = mapInPlace ArithmeticOperations.intNotQ clContext workGroupSize
+
+        let prefixSum = PrefixSum.standardExcludeInPlace clContext workGroupSize
+
+        let scatter = Scatter.lastOccurrence clContext workGroupSize
+
+        fun (queue: MailboxProcessor<_>) allocationMode (excludeBitmap: ClArray<int>) (inputArray: ClArray<'a>) ->
+
+            invert queue excludeBitmap
+
+            let length = (prefixSum queue excludeBitmap).ToHostAndFree queue
+
+            if length = 0 then
+                None
+            else
+                let result = clContext.CreateClArrayWithSpecificAllocationMode(allocationMode, length)
+
+                scatter queue excludeBitmap inputArray result
+
+                Some result

src/GraphBLAS-sharp.Backend/Common/Map.fs

@@ -40,6 +40,36 @@ module Map =
 
             result
 
+    /// <summary>
+    /// Changes elements of the input array, applying the given function
+    /// to each element of the array.
+    /// </summary>
+    /// <param name="op">The function to transform elements of the array.</param>
+    /// <param name="clContext">OpenCL context.</param>
+    /// <param name="workGroupSize">Should be a power of 2 and greater than 1.</param>
+    let mapInPlace<'a> (op: Expr<'a -> 'a>) (clContext: ClContext) workGroupSize =
+
+        let map =
+            <@ fun (ndRange: Range1D) lenght (inputArray: ClArray<'a>) ->
+
+                let gid = ndRange.GlobalID0
+
+                if gid < lenght then
+                    inputArray.[gid] <- (%op) inputArray.[gid] @>
+
+        let kernel = clContext.Compile map
+
+        fun (processor: MailboxProcessor<_>) (inputArray: ClArray<'a>) ->
+
+            let ndRange =
+                Range1D.CreateValid(inputArray.Length, workGroupSize)
+
+            let kernel = kernel.GetKernel()
+
+            processor.Post(Msg.MsgSetArguments(fun () -> kernel.KernelFunc ndRange inputArray.Length inputArray))
+
+            processor.Post(Msg.CreateRunMsg<_, _>(kernel))
+
     /// <summary>
     /// Builds a new array whose elements are the results of applying the given function
     /// to the corresponding pairs of values, where the first element of pair is from the given array

src/GraphBLAS-sharp.Backend/GraphBLAS-sharp.Backend.fsproj

@@ -69,6 +69,7 @@
     <Compile Include="Operations/Operations.fs" />
 
     <Compile Include="Algorithms/BFS.fs" />
+    <Compile Include="Algorithms/MSBFS.fs" />
     <Compile Include="Algorithms/SSSP.fs" />
     <Compile Include="Algorithms/PageRank.fs" />
     <Compile Include="Algorithms/Algorithms.fs" />

src/GraphBLAS-sharp.Backend/Matrix/COO/Intersect.fs

@@ -1,13 +1,12 @@
 namespace GraphBLAS.FSharp.Backend.Matrix.COO
 
 open Brahma.FSharp
-open GraphBLAS.FSharp
 open GraphBLAS.FSharp.Objects
 open GraphBLAS.FSharp.Backend.Quotes
 open GraphBLAS.FSharp.Objects.ClContextExtensions
 
 module internal Intersect =
-    let findIntersectionByKeys (clContext: ClContext) workGroupSize =
+    let findKeysIntersection (clContext: ClContext) workGroupSize =
 
         let findIntersection =
             <@ fun (ndRange: Range1D) (leftNNZ: int) (rightNNZ: int) (leftRows: ClArray<int>) (leftColumns: ClArray<int>) (rightRows: ClArray<int>) (rightColumns: ClArray<int>) (bitmap: ClArray<int>) ->
@@ -17,15 +16,9 @@ module internal Intersect =
 
                 if gid < bitmapSize then
 
-                    let row, column, rows, columns =
-                        if rightNNZ >= leftNNZ then
-                            leftRows.[gid], leftColumns.[gid], leftRows, leftColumns
-                        else
-                            rightRows.[gid], rightColumns.[gid], rightRows, rightColumns
+                    let index: uint64 = ((uint64 leftRows.[gid]) <<< 32) ||| (uint64 leftColumns.[gid])
 
-                    let index: uint64 = ((uint64 row) <<< 32) ||| (uint64 column)
-
-                    let intersect = (%Search.Bin.existsByKey2D) bitmapSize index rows columns
+                    let intersect = (%Search.Bin.existsByKey2D) bitmapSize index rightRows rightColumns
 
                     if intersect then
                         bitmap.[gid] <- 1
@@ -37,7 +30,7 @@ module internal Intersect =
 
         fun (processor: MailboxProcessor<_>) allocationMode (leftMatrix: ClMatrix.COO<'a>) (rightMatrix: ClMatrix.COO<'b>) ->
 
-            let bitmapSize = min leftMatrix.NNZ rightMatrix.NNZ
+            let bitmapSize = leftMatrix.NNZ
 
             let bitmap = clContext.CreateClArrayWithSpecificAllocationMode(allocationMode, bitmapSize)
 

src/GraphBLAS-sharp.Backend/Matrix/COO/Map.fs

@@ -1,6 +1,5 @@
 namespace GraphBLAS.FSharp.Backend.Matrix.COO
 
-open System
 open Brahma.FSharp
 open GraphBLAS.FSharp.Backend.Matrix
 open GraphBLAS.FSharp.Backend.Quotes

src/GraphBLAS-sharp.Backend/Matrix/COO/Matrix.fs

@@ -8,6 +8,7 @@ open GraphBLAS.FSharp.Objects
 open GraphBLAS.FSharp.Objects.ClMatrix
 open GraphBLAS.FSharp.Objects.ClCellExtensions
 open GraphBLAS.FSharp.Objects.ArraysExtensions
+open GraphBLAS.FSharp.Objects.ClContextExtensions
 
 module Matrix =
     /// <summary>
@@ -218,10 +219,35 @@ module Matrix =
             |> transposeInPlace queue
 
     /// <summary>
-    /// Build a bitmap. Maps non-zero elements of the matrix with minimum nnz to 1
-    /// if the second matrix has non zero element under the same row and column pair, otherwise 0.
+    /// Builds a bitmap. Maps non-zero elements of the left matrix
+    /// to 1 if the right matrix has non zero element under the same row and column pair, otherwise 0.
     /// </summary>
     /// <param name="clContext">OpenCL context.</param>
     /// <param name="workGroupSize">Should be a power of 2 and greater than 1.</param>
     let findIntersectionByKeys (clContext: ClContext) workGroupSize =
-        Intersect.findIntersectionByKeys clContext workGroupSize
+        Intersect.findKeysIntersection clContext workGroupSize
+
+    /// <summary>
+    /// Merges two disjoint matrices of the same size.
+    /// </summary>
+    /// <remarks>
+    /// Matrices should have the same number of rows and columns. <br/>
+    /// Matrices should not have non zero values with the same index.
+    /// </remarks>
+    /// <param name="clContext">OpenCL context.</param>
+    /// <param name="workGroupSize">Should be a power of 2 and greater than 1.</param>
+    let mergeDisjoint (clContext: ClContext) workGroupSize = Merge.runDisjoint clContext workGroupSize
+
+    let ofList (clContext: ClContext) allocationMode rowCount columnCount (elements: (int * int * 'a) list) =
+        let rows, columns, values =
+            elements
+            |> Array.ofList
+            |> Array.sortBy (fun (x, _, _) -> x)
+            |> Array.unzip3
+
+        { Context = clContext
+          Rows = clContext.CreateClArrayWithSpecificAllocationMode(allocationMode, rows)
+          Columns = clContext.CreateClArrayWithSpecificAllocationMode(allocationMode, columns)
+          Values = clContext.CreateClArrayWithSpecificAllocationMode(allocationMode, values)
+          RowCount = rowCount
+          ColumnCount = columnCount }