diff --git a/QuadTree.Tests/Tests.Matrix.fs b/QuadTree.Tests/Tests.Matrix.fs
index 816b4a2..74e4c42 100644
--- a/QuadTree.Tests/Tests.Matrix.fs
+++ b/QuadTree.Tests/Tests.Matrix.fs
@@ -642,3 +642,1686 @@ let ``Fold sum`` () =
     let actual = foldAssociative op_add None m1 |> Option.get
 
     Assert.Equal(expected, actual)
+
+[<Fact>]
+let ``fromCoordinateList with out-of-range coordinates throws exception`` () =
+    let coo =
+        CoordinateList(6UL<nrows>, 6UL<ncols>, [ (9UL<rowindex>, 9UL<colindex>, 13) ])
+
+    Assert.Throws<Exception>(fun () -> fromCoordinateList coo |> ignore)
+
+[<Fact>]
+let ``fromCoordinateList with unsorted coordinates works correctly`` () =
+    let coo =
+        CoordinateList(
+            7UL<nrows>,
+            7UL<ncols>,
+            [ (6UL<rowindex>, 6UL<colindex>, 10)
+              (1UL<rowindex>, 2UL<colindex>, 8)
+              (1UL<rowindex>, 1UL<colindex>, 100) ]
+        )
+
+    let matrix = fromCoordinateList coo
+    let result = toCoordinateList matrix
+
+    Assert.Equal(
+        CoordinateList(
+            7UL<nrows>,
+            7UL<ncols>,
+            [ (1UL<rowindex>, 1UL<colindex>, 100)
+              (1UL<rowindex>, 2UL<colindex>, 8)
+              (6UL<rowindex>, 6UL<colindex>, 10) ]
+        ),
+        result
+    )
+
+[<Fact>]
+let ``fromCoordinateList with duplicate indices returns the last of them`` () =
+    let coo =
+        CoordinateList(
+            3UL<nrows>,
+            3UL<ncols>,
+            [ (1UL<rowindex>, 1UL<colindex>, 33); (1UL<rowindex>, 1UL<colindex>, 100) ]
+        )
+
+    let matrix = fromCoordinateList coo
+    let result = toCoordinateList matrix
+    Assert.Equal(CoordinateList(3UL<nrows>, 3UL<ncols>, [ (1UL<rowindex>, 1UL<colindex>, 100) ]), result)
+
+[<Fact>]
+let ``fromCoordinateList with zero size throws Exception`` () =
+    let coo =
+        CoordinateList(
+            0UL<nrows>,
+            0UL<ncols>,
+            [ (33UL<rowindex>, 33UL<colindex>, 33); (39UL<rowindex>, 39UL<colindex>, 1) ]
+        )
+
+    Assert.Throws<Exception>(fun () -> fromCoordinateList coo |> ignore)
+
+[<Fact>]
+let ``map works on square matrix`` () =
+    let matrix =
+        fromCoordinateList (
+            CoordinateList(
+                4UL<nrows>,
+                4UL<ncols>,
+                [ (1UL<rowindex>, 1UL<colindex>, 22); (2UL<rowindex>, 3UL<colindex>, 37) ]
+            )
+        )
+
+    let result =
+        map matrix (fun (x: int option) ->
+            match x with
+            | Some(v) -> Some(v + 9)
+            | None -> None)
+
+    let coo = toCoordinateList result
+
+    Assert.Equal(
+        CoordinateList(
+            4UL<nrows>,
+            4UL<ncols>,
+            [ (1UL<rowindex>, 1UL<colindex>, 31); (2UL<rowindex>, 3UL<colindex>, 46) ]
+        ),
+        coo
+    )
+
+[<Fact>]
+let ``map works on rectangular matrix`` () =
+    let matrix =
+        fromCoordinateList (
+            CoordinateList(
+                5UL<nrows>,
+                6UL<ncols>,
+                [ (1UL<rowindex>, 1UL<colindex>, 21); (4UL<rowindex>, 3UL<colindex>, 36) ]
+            )
+        )
+
+    let result =
+        map matrix (fun (x: int option) ->
+            match x with
+            | Some(v) -> Some(v / 3)
+            | None -> None)
+
+    let coo = toCoordinateList result
+
+    Assert.Equal(
+        CoordinateList(
+            5UL<nrows>,
+            6UL<ncols>,
+            [ (1UL<rowindex>, 1UL<colindex>, 7); (4UL<rowindex>, 3UL<colindex>, 12) ]
+        ),
+        coo
+    )
+
+[<Fact>]
+let ``map on empty matrix returns empty matrix`` () =
+    let matrix = fromCoordinateList (CoordinateList(0UL<nrows>, 0UL<ncols>, []))
+
+    let result =
+        map matrix (fun (x: int option) ->
+            match x with
+            | Some(v) -> Some(v * 5)
+            | None -> None)
+
+    let coo = toCoordinateList result
+    Assert.Equal(CoordinateList(0UL<nrows>, 0UL<ncols>, []), coo)
+
+[<Fact>]
+let ``map with function that turns all the elements into zeros`` () =
+    let matrix =
+        fromCoordinateList (
+            CoordinateList(
+                5UL<nrows>,
+                7UL<ncols>,
+                [ (1UL<rowindex>, 1UL<colindex>, 5)
+                  (4UL<rowindex>, 1UL<colindex>, 17)
+                  (4UL<rowindex>, 6UL<colindex>, 33) ]
+            )
+        )
+
+    let result = map matrix (fun _ -> Some(0))
+    let coo = toCoordinateList result
+
+    let expectedData =
+        [ for i in 0UL .. 4UL do
+              for j in 0UL .. 6UL -> (i * 1UL<rowindex>, j * 1UL<colindex>, 0) ]
+
+    let expected = CoordinateList(5UL<nrows>, 7UL<ncols>, expectedData)
+    Assert.Equal(expected, coo)
+
+[<Fact>]
+let ``map with function that turns all the elements to None`` () =
+    let matrix =
+        fromCoordinateList (
+            CoordinateList(
+                5UL<nrows>,
+                7UL<ncols>,
+                [ (1UL<rowindex>, 1UL<colindex>, 5)
+                  (4UL<rowindex>, 1UL<colindex>, 17)
+                  (4UL<rowindex>, 6UL<colindex>, 33) ]
+            )
+        )
+
+    let result = map matrix (fun _ -> None)
+    let coo = toCoordinateList result
+    Assert.Equal(CoordinateList(5UL<nrows>, 7UL<ncols>, []), coo)
+
+[<Fact>]
+let ``map can change type from int to string`` () =
+    let matrix =
+        fromCoordinateList (
+            CoordinateList(
+                3UL<nrows>,
+                5UL<ncols>,
+                [ (2UL<rowindex>, 1UL<colindex>, 17); (2UL<rowindex>, 4UL<colindex>, 33) ]
+            )
+        )
+
+    let result =
+        map matrix (fun (x: int option) ->
+            match x with
+            | Some v -> Some(sprintf "str %d" v)
+            | None -> None)
+
+    let coo = toCoordinateList result
+
+    Assert.Equal(
+        CoordinateList(
+            3UL<nrows>,
+            5UL<ncols>,
+            [ (2UL<rowindex>, 1UL<colindex>, "str 17")
+              (2UL<rowindex>, 4UL<colindex>, "str 33") ]
+        ),
+        coo
+    )
+
+[<Fact>]
+let ``mapi works with row index and col index on square matrix`` () =
+    let coo =
+        CoordinateList(
+            4UL<nrows>,
+            4UL<ncols>,
+            [ (1UL<rowindex>, 1UL<colindex>, 12); (2UL<rowindex>, 3UL<colindex>, 13) ]
+        )
+
+    let matrix = fromCoordinateList coo
+
+    let result =
+        mapi matrix (fun (i: uint64<rowindex>) (j: uint64<colindex>) (x: int option) ->
+            match x with
+            | Some v -> Some(v * int (i) + int (j))
+            | None -> None)
+
+    let actual = toCoordinateList result
+
+    Assert.Equal(
+        CoordinateList(
+            4UL<nrows>,
+            4UL<ncols>,
+            [ (1UL<rowindex>, 1UL<colindex>, 13); (2UL<rowindex>, 3UL<colindex>, 29) ]
+        ),
+        actual
+    )
+
+[<Fact>]
+let ``mapi works with row index and col index on rectangular matrix`` () =
+    let coo =
+        CoordinateList(
+            3UL<nrows>,
+            5UL<ncols>,
+            [ (1UL<rowindex>, 2UL<colindex>, 15); (2UL<rowindex>, 2UL<colindex>, 13) ]
+        )
+
+    let matrix = fromCoordinateList coo
+
+    let result =
+        mapi matrix (fun (i: uint64<rowindex>) (j: uint64<colindex>) (x: int option) ->
+            match x with
+            | Some v -> Some(v * int (j) + int (i))
+            | None -> None)
+
+    let actual = toCoordinateList result
+
+    Assert.Equal(
+        CoordinateList(
+            3UL<nrows>,
+            5UL<ncols>,
+            [ (1UL<rowindex>, 2UL<colindex>, 31); (2UL<rowindex>, 2UL<colindex>, 28) ]
+        ),
+        actual
+    )
+
+[<Fact>]
+let ``mapi on empty matrix returns empty matrix`` () =
+    let matrix = fromCoordinateList (CoordinateList(0UL<nrows>, 0UL<ncols>, []))
+
+    let result =
+        mapi matrix (fun (i: uint64<rowindex>) (j: uint64<colindex>) (x: int option) ->
+            match x with
+            | Some v -> Some((v + int (1)) * 3 + 2 * int (j))
+            | None -> None)
+
+    let actual = toCoordinateList result
+    Assert.Equal(CoordinateList(0UL<nrows>, 0UL<ncols>, []), actual)
+
+[<Fact>]
+let ``mapi with special function returns empty matrix`` () =
+    let coo =
+        CoordinateList(
+            5UL<nrows>,
+            7UL<ncols>,
+            [ (1UL<rowindex>, 2UL<colindex>, 15); (4UL<rowindex>, 5UL<colindex>, 13) ]
+        )
+
+    let matrix = fromCoordinateList coo
+
+    let result =
+        mapi matrix (fun (i: uint64<rowindex>) (j: uint64<colindex>) (x: int option) ->
+            match x with
+            | Some v ->
+                match (int (i) + int (j)) % 2 with
+                | 0 -> Some(v * 2)
+                | _ -> None
+            | None -> None)
+
+    let actual = toCoordinateList result
+    Assert.Equal(CoordinateList(5UL<nrows>, 7UL<ncols>, []), actual)
+
+[<Fact>]
+let ``mapi works with function does not depend on the indexes`` () =
+    let coo =
+        CoordinateList(
+            5UL<nrows>,
+            7UL<ncols>,
+            [ (1UL<rowindex>, 2UL<colindex>, 15); (4UL<rowindex>, 5UL<colindex>, 13) ]
+        )
+
+    let matrix = fromCoordinateList coo
+
+    let result =
+        mapi matrix (fun (i: uint64<rowindex>) (j: uint64<colindex>) (x: int option) ->
+            match x with
+            | Some v -> Some(v * 2)
+            | None -> None)
+
+    let actual = toCoordinateList result
+
+    Assert.Equal(
+        CoordinateList(
+            5UL<nrows>,
+            7UL<ncols>,
+            [ (1UL<rowindex>, 2UL<colindex>, 30); (4UL<rowindex>, 5UL<colindex>, 26) ]
+        ),
+        actual
+    )
+
+[<Fact>]
+let ``slice returns error when row start is negative`` () =
+    let matrix =
+        fromCoordinateList (
+            CoordinateList(
+                5UL<nrows>,
+                6UL<ncols>,
+                [ (2UL<rowindex>, 2UL<colindex>, 3); (3UL<rowindex>, 4UL<colindex>, 17) ]
+            )
+        )
+
+    match slice matrix -1 4 2 3 with
+    | Result.Ok _ -> Assert.Fail("Expected Error")
+    | Result.Error msg -> Assert.Equal("Start row should be >= 0", msg)
+
+[<Fact>]
+let ``slice returns error when row end is negative`` () =
+    let matrix =
+        fromCoordinateList (
+            CoordinateList(
+                5UL<nrows>,
+                6UL<ncols>,
+                [ (2UL<rowindex>, 2UL<colindex>, 3); (3UL<rowindex>, 4UL<colindex>, 17) ]
+            )
+        )
+
+    match slice matrix 1 -4 2 3 with
+    | Result.Ok _ -> Assert.Fail("Expected Error")
+    | Result.Error msg -> Assert.Equal("End row should be >= 0", msg)
+
+[<Fact>]
+let ``slice returns error when col start is negative`` () =
+    let matrix =
+        fromCoordinateList (
+            CoordinateList(
+                5UL<nrows>,
+                6UL<ncols>,
+                [ (2UL<rowindex>, 2UL<colindex>, 3); (3UL<rowindex>, 4UL<colindex>, 17) ]
+            )
+        )
+
+    match slice matrix 1 4 -2 3 with
+    | Result.Ok _ -> Assert.Fail("Expected Error")
+    | Result.Error msg -> Assert.Equal("Start column should be >= 0", msg)
+
+[<Fact>]
+let ``slice returns error when col end is negative`` () =
+    let matrix =
+        fromCoordinateList (
+            CoordinateList(
+                5UL<nrows>,
+                6UL<ncols>,
+                [ (2UL<rowindex>, 2UL<colindex>, 3); (3UL<rowindex>, 4UL<colindex>, 17) ]
+            )
+        )
+
+    match slice matrix 1 4 2 -3 with
+    | Result.Ok _ -> Assert.Fail("Expected Error")
+    | Result.Error msg -> Assert.Equal("End column should be >= 0", msg)
+
+[<Fact>]
+let ``slice returns error when row start is out of range`` () =
+    let matrix =
+        fromCoordinateList (
+            CoordinateList(
+                5UL<nrows>,
+                6UL<ncols>,
+                [ (2UL<rowindex>, 2UL<colindex>, 3); (3UL<rowindex>, 4UL<colindex>, 17) ]
+            )
+        )
+
+    match slice matrix 6 4 2 3 with
+    | Result.Ok _ -> Assert.Fail("Expected Error")
+    | Result.Error msg -> Assert.Equal("Start row is out of matrix length", msg)
+
+[<Fact>]
+let ``slice returns error when row end is out of range`` () =
+    let matrix =
+        fromCoordinateList (
+            CoordinateList(
+                5UL<nrows>,
+                6UL<ncols>,
+                [ (2UL<rowindex>, 2UL<colindex>, 3); (3UL<rowindex>, 4UL<colindex>, 17) ]
+            )
+        )
+
+    match slice matrix 1 10 2 3 with
+    | Result.Ok _ -> Assert.Fail("Expected Error")
+    | Result.Error msg -> Assert.Equal("End row is out of matrix length", msg)
+
+[<Fact>]
+let ``slice returns error when col start is out of range`` () =
+    let matrix =
+        fromCoordinateList (
+            CoordinateList(
+                5UL<nrows>,
+                6UL<ncols>,
+                [ (2UL<rowindex>, 2UL<colindex>, 3); (3UL<rowindex>, 4UL<colindex>, 17) ]
+            )
+        )
+
+    match slice matrix 1 4 10 3 with
+    | Result.Ok _ -> Assert.Fail("Expected Error")
+    | Result.Error msg -> Assert.Equal("Start column is out of matrix length", msg)
+
+[<Fact>]
+let ``slice returns error when col end is out of range`` () =
+    let matrix =
+        fromCoordinateList (
+            CoordinateList(
+                5UL<nrows>,
+                6UL<ncols>,
+                [ (2UL<rowindex>, 2UL<colindex>, 3); (3UL<rowindex>, 4UL<colindex>, 17) ]
+            )
+        )
+
+    match slice matrix 1 2 2 10 with
+    | Result.Ok _ -> Assert.Fail("Expected Error")
+    | Result.Error msg -> Assert.Equal("End column is out of matrix length", msg)
+
+[<Fact>]
+let ``slice returns error when row end is less than row start`` () =
+    let matrix =
+        fromCoordinateList (
+            CoordinateList(
+                5UL<nrows>,
+                6UL<ncols>,
+                [ (2UL<rowindex>, 2UL<colindex>, 3); (3UL<rowindex>, 4UL<colindex>, 17) ]
+            )
+        )
+
+    match slice matrix 2 1 2 3 with
+    | Result.Ok _ -> Assert.Fail("Expected Error")
+    | Result.Error msg -> Assert.Equal("Start row should be <= end row", msg)
+
+[<Fact>]
+let ``slice returns error when col end is less than col start`` () =
+    let matrix =
+        fromCoordinateList (
+            CoordinateList(
+                5UL<nrows>,
+                6UL<ncols>,
+                [ (2UL<rowindex>, 2UL<colindex>, 3); (3UL<rowindex>, 4UL<colindex>, 17) ]
+            )
+        )
+
+    match slice matrix 1 2 3 2 with
+    | Result.Ok _ -> Assert.Fail("Expected Error")
+    | Result.Error msg -> Assert.Equal("Start column should be <= end column", msg)
+
+[<Fact>]
+let ``slice returns correct square submatrix`` () =
+    let matrix =
+        fromCoordinateList (
+            CoordinateList(
+                7UL<nrows>,
+                7UL<ncols>,
+                [ (2UL<rowindex>, 2UL<colindex>, 33); (5UL<rowindex>, 5UL<colindex>, 28) ]
+            )
+        )
+
+    match slice matrix 1 3 1 3 with
+    | Result.Ok result ->
+        let coo = toCoordinateList result
+        Assert.Equal(CoordinateList(3UL<nrows>, 3UL<ncols>, [ (1UL<rowindex>, 1UL<colindex>, 33) ]), coo)
+    | Result.Error msg -> Assert.Fail(msg)
+
+[<Fact>]
+let ``slice returns correct rectangular submatrix`` () =
+    let matrix =
+        fromCoordinateList (
+            CoordinateList(
+                5UL<nrows>,
+                9UL<ncols>,
+                [ (3UL<rowindex>, 7UL<colindex>, 33); (1UL<rowindex>, 2UL<colindex>, 28) ]
+            )
+        )
+
+    match slice matrix 2 4 5 8 with
+    | Result.Ok result ->
+        let coo = toCoordinateList result
+        Assert.Equal(CoordinateList(3UL<nrows>, 4UL<ncols>, [ (1UL<rowindex>, 2UL<colindex>, 33) ]), coo)
+    | Result.Error msg -> Assert.Fail(msg)
+
+[<Fact>]
+let ``slice returns empty matrix`` () =
+    let matrix =
+        fromCoordinateList (
+            CoordinateList(
+                7UL<nrows>,
+                9UL<ncols>,
+                [ (3UL<rowindex>, 3UL<colindex>, 33); (1UL<rowindex>, 2UL<colindex>, 28) ]
+            )
+        )
+
+    match slice matrix 4 6 5 8 with
+    | Result.Ok result ->
+        let coo = toCoordinateList result
+        Assert.Equal(CoordinateList(3UL<nrows>, 4UL<ncols>, []), coo)
+    | Result.Error msg -> Assert.Fail(msg)
+
+[<Fact>]
+let ``slice returns single submatrix`` () =
+    let matrix =
+        fromCoordinateList (
+            CoordinateList(
+                7UL<nrows>,
+                9UL<ncols>,
+                [ (3UL<rowindex>, 3UL<colindex>, 33); (1UL<rowindex>, 2UL<colindex>, 28) ]
+            )
+        )
+
+    match slice matrix 1 1 2 2 with
+    | Result.Ok result ->
+        let coo = toCoordinateList result
+        Assert.Equal(CoordinateList(1UL<nrows>, 1UL<ncols>, [ (0UL<rowindex>, 0UL<colindex>, 28) ]), coo)
+    | Result.Error msg -> Assert.Fail(msg)
+
+[<Fact>]
+let ``slice returns correct submatrix equals to matrix`` () =
+    let matrix =
+        fromCoordinateList (
+            CoordinateList(
+                7UL<nrows>,
+                9UL<ncols>,
+                [ (3UL<rowindex>, 3UL<colindex>, 33); (1UL<rowindex>, 2UL<colindex>, 28) ]
+            )
+        )
+
+    match slice matrix 0 6 0 8 with
+    | Result.Ok result ->
+        let coo = toCoordinateList result
+
+        Assert.Equal(
+            CoordinateList(
+                7UL<nrows>,
+                9UL<ncols>,
+                [ (1UL<rowindex>, 2UL<colindex>, 28); (3UL<rowindex>, 3UL<colindex>, 33) ]
+            ),
+            coo
+        )
+    | Result.Error msg -> Assert.Fail(msg)
+
+[<Fact>]
+let ``slice returns correct submatrix when row start of submatrix equals to row start of matrix`` () =
+    let matrix =
+        fromCoordinateList (
+            CoordinateList(
+                7UL<nrows>,
+                7UL<ncols>,
+                [ (0UL<rowindex>, 0UL<colindex>, 10)
+                  (3UL<rowindex>, 3UL<colindex>, 33)
+                  (6UL<rowindex>, 6UL<colindex>, 6) ]
+            )
+        )
+
+    match slice matrix 0 4 0 4 with
+    | Result.Ok result ->
+        let coo = toCoordinateList result
+
+        Assert.Equal(
+            CoordinateList(
+                5UL<nrows>,
+                5UL<ncols>,
+                [ (0UL<rowindex>, 0UL<colindex>, 10); (3UL<rowindex>, 3UL<colindex>, 33) ]
+            ),
+            coo
+        )
+    | Result.Error msg -> Assert.Fail(msg)
+
+[<Fact>]
+let ``slice returns correct submatrix when row end of submatrix equals to row end of matrix`` () =
+    let matrix =
+        fromCoordinateList (
+            CoordinateList(
+                7UL<nrows>,
+                7UL<ncols>,
+                [ (0UL<rowindex>, 0UL<colindex>, 10)
+                  (3UL<rowindex>, 3UL<colindex>, 33)
+                  (6UL<rowindex>, 6UL<colindex>, 6) ]
+            )
+        )
+
+    match slice matrix 2 6 2 4 with
+    | Result.Ok result ->
+        let coo = toCoordinateList result
+        Assert.Equal(CoordinateList(5UL<nrows>, 3UL<ncols>, [ (1UL<rowindex>, 1UL<colindex>, 33) ]), coo)
+    | Result.Error msg -> Assert.Fail(msg)
+
+[<Fact>]
+let ``slice returns correct submatrix when col start of submatrix equals to col start of matrix`` () =
+    let matrix =
+        fromCoordinateList (
+            CoordinateList(
+                7UL<nrows>,
+                7UL<ncols>,
+                [ (0UL<rowindex>, 0UL<colindex>, 10)
+                  (3UL<rowindex>, 3UL<colindex>, 33)
+                  (6UL<rowindex>, 6UL<colindex>, 6) ]
+            )
+        )
+
+    match slice matrix 2 5 0 6 with
+    | Result.Ok result ->
+        let coo = toCoordinateList result
+        Assert.Equal(CoordinateList(4UL<nrows>, 7UL<ncols>, [ (1UL<rowindex>, 3UL<colindex>, 33) ]), coo)
+    | Result.Error msg -> Assert.Fail(msg)
+
+[<Fact>]
+let ``slice returns correct submatrix when col end of submatrix equals to col end of matrix`` () =
+    let matrix =
+        fromCoordinateList (
+            CoordinateList(
+                7UL<nrows>,
+                7UL<ncols>,
+                [ (0UL<rowindex>, 0UL<colindex>, 10)
+                  (3UL<rowindex>, 3UL<colindex>, 33)
+                  (6UL<rowindex>, 6UL<colindex>, 6) ]
+            )
+        )
+
+    match slice matrix 2 4 2 6 with
+    | Result.Ok result ->
+        let coo = toCoordinateList result
+        Assert.Equal(CoordinateList(3UL<nrows>, 5UL<ncols>, [ (1UL<rowindex>, 1UL<colindex>, 33) ]), coo)
+    | Result.Error msg -> Assert.Fail(msg)
+
+[<Fact>]
+let ``slice returns single column`` () =
+    let matrix =
+        fromCoordinateList (
+            CoordinateList(
+                10UL<nrows>,
+                10UL<ncols>,
+                [ (1UL<rowindex>, 3UL<colindex>, 1)
+                  (2UL<rowindex>, 2UL<colindex>, 2)
+                  (5UL<rowindex>, 7UL<colindex>, 3) ]
+            )
+        )
+
+    match slice matrix 1 6 2 2 with
+    | Result.Ok result ->
+        let coo = toCoordinateList result
+        Assert.Equal(CoordinateList(6UL<nrows>, 1UL<ncols>, [ (1UL<rowindex>, 0UL<colindex>, 2) ]), coo)
+    | Result.Error msg -> Assert.Fail(msg)
+
+[<Fact>]
+let ``slice returns single row`` () =
+    let matrix =
+        fromCoordinateList (
+            CoordinateList(
+                10UL<nrows>,
+                10UL<ncols>,
+                [ (1UL<rowindex>, 3UL<colindex>, 1)
+                  (2UL<rowindex>, 2UL<colindex>, 2)
+                  (5UL<rowindex>, 7UL<colindex>, 3) ]
+            )
+        )
+
+    match slice matrix 5 5 3 9 with
+    | Result.Ok result ->
+        let coo = toCoordinateList result
+        Assert.Equal(CoordinateList(1UL<nrows>, 7UL<ncols>, [ (0UL<rowindex>, 4UL<colindex>, 3) ]), coo)
+    | Result.Error msg -> Assert.Fail(msg)
+
+[<Fact>]
+let ``let reduceRows sum on square power of two matrix`` () =
+    let coo =
+        CoordinateList(
+            2UL<nrows>,
+            2UL<ncols>,
+            [ (0UL<rowindex>, 0UL<colindex>, 12)
+              (0UL<rowindex>, 1UL<colindex>, 15)
+              (1UL<rowindex>, 0UL<colindex>, 17)
+              (1UL<rowindex>, 1UL<colindex>, 3) ]
+        )
+
+    let matrix = fromCoordinateList coo
+
+    let add x y =
+        match x, y with
+        | None, None -> None
+        | None, Some(b) -> Some(b)
+        | Some(a), None -> Some(a)
+        | Some(a), Some(b) -> Some(a + b)
+
+    let result = reduceRows add matrix
+    let vectorCoordinates = Vector.toCoordinateList result
+
+    let expected =
+        Vector.CoordinateList(2UL<Vector.dataLength>, [ (0UL<Vector.index>, 27); (1UL<Vector.index>, 20) ])
+
+    Assert.Equal(expected, vectorCoordinates)
+
+[<Fact>]
+let ``let reduceRows sum on square power of two matrix with empty row`` () =
+    let coo =
+        CoordinateList(
+            2UL<nrows>,
+            2UL<ncols>,
+            [ (1UL<rowindex>, 0UL<colindex>, 17); (1UL<rowindex>, 1UL<colindex>, 3) ]
+        )
+
+    let matrix = fromCoordinateList coo
+
+    let add x y =
+        match x, y with
+        | None, None -> None
+        | None, Some(b) -> Some(b)
+        | Some(a), None -> Some(a)
+        | Some(a), Some(b) -> Some(a + b)
+
+    let result = reduceRows add matrix
+    let vectorCoordinates = Vector.toCoordinateList result
+
+    let expected =
+        Vector.CoordinateList(2UL<Vector.dataLength>, [ (1UL<Vector.index>, 20) ])
+
+    Assert.Equal(expected, vectorCoordinates)
+
+[<Fact>]
+let ``let reduceRows sum on square not power of two matrix`` () =
+    let coo =
+        CoordinateList(
+            3UL<nrows>,
+            3UL<ncols>,
+            [ (0UL<rowindex>, 0UL<colindex>, 5)
+              (0UL<rowindex>, 1UL<colindex>, 7)
+              (0UL<rowindex>, 2UL<colindex>, 9)
+              (1UL<rowindex>, 1UL<colindex>, 11)
+              (1UL<rowindex>, 2UL<colindex>, 13)
+              (2UL<rowindex>, 0UL<colindex>, 15)
+              (2UL<rowindex>, 1UL<colindex>, 17) ]
+        )
+
+    let matrix = fromCoordinateList coo
+
+    let add x y =
+        match x, y with
+        | None, None -> None
+        | None, Some(b) -> Some(b)
+        | Some(a), None -> Some(a)
+        | Some(a), Some(b) -> Some(a + b)
+
+    let result = reduceRows add matrix
+    let vectorCoordinates = Vector.toCoordinateList result
+
+    let expected =
+        Vector.CoordinateList(
+            3UL<Vector.dataLength>,
+            [ (0UL<Vector.index>, 21); (1UL<Vector.index>, 24); (2UL<Vector.index>, 32) ]
+        )
+
+    Assert.Equal(expected, vectorCoordinates)
+
+[<Fact>]
+let ``let reduceRows mul on square not power of two matrix`` () =
+    let coo =
+        CoordinateList(
+            3UL<nrows>,
+            3UL<ncols>,
+            [ (0UL<rowindex>, 0UL<colindex>, 5)
+              (0UL<rowindex>, 1UL<colindex>, 7)
+              (0UL<rowindex>, 2UL<colindex>, 9)
+              (1UL<rowindex>, 1UL<colindex>, 11)
+              (1UL<rowindex>, 2UL<colindex>, 13)
+              (2UL<rowindex>, 0UL<colindex>, 15)
+              (2UL<rowindex>, 1UL<colindex>, 17) ]
+        )
+
+    let matrix = fromCoordinateList coo
+
+    let mul x y =
+        match x, y with
+        | None, None -> None
+        | None, Some(b) -> Some(b)
+        | Some(a), None -> Some(a)
+        | Some(a), Some(b) -> Some(a * b)
+
+    let result = reduceRows mul matrix
+    let vectorCoordinates = Vector.toCoordinateList result
+
+    let expected =
+        Vector.CoordinateList(
+            3UL<Vector.dataLength>,
+            [ (0UL<Vector.index>, 315); (1UL<Vector.index>, 143); (2UL<Vector.index>, 255) ]
+        )
+
+    Assert.Equal(expected, vectorCoordinates)
+
+[<Fact>]
+let ``let reduceRows sum on rectangular matrix`` () =
+    let coo =
+        CoordinateList(
+            2UL<nrows>,
+            3UL<ncols>,
+            [ (0UL<rowindex>, 0UL<colindex>, 5)
+              (0UL<rowindex>, 1UL<colindex>, 7)
+              (0UL<rowindex>, 2UL<colindex>, 9)
+              (1UL<rowindex>, 1UL<colindex>, 11)
+              (1UL<rowindex>, 2UL<colindex>, 13) ]
+        )
+
+    let matrix = fromCoordinateList coo
+
+    let sum x y =
+        match x, y with
+        | None, None -> None
+        | None, Some(b) -> Some(b)
+        | Some(a), None -> Some(a)
+        | Some(a), Some(b) -> Some(a + b)
+
+    let result = reduceRows sum matrix
+    let vectorCoordinates = Vector.toCoordinateList result
+
+    let expected =
+        Vector.CoordinateList(2UL<Vector.dataLength>, [ (0UL<Vector.index>, 21); (1UL<Vector.index>, 24) ])
+
+    Assert.Equal(expected, vectorCoordinates)
+
+[<Fact>]
+let ``let reduceRows sum on empty matrix`` () =
+    let coo = CoordinateList(2UL<nrows>, 3UL<ncols>, [])
+    let matrix = fromCoordinateList coo
+
+    let sum x y =
+        match x, y with
+        | None, None -> None
+        | None, Some(b) -> Some(b)
+        | Some(a), None -> Some(a)
+        | Some(a), Some(b) -> Some(a + b)
+
+    let result = reduceRows sum matrix
+    let vectorCoordinates = Vector.toCoordinateList result
+    let expected = Vector.CoordinateList(2UL<Vector.dataLength>, [])
+    Assert.Equal(expected, vectorCoordinates)
+
+[<Fact>]
+let ``let reduceRows mul on single matrix`` () =
+    let coo =
+        CoordinateList(1UL<nrows>, 1UL<ncols>, [ (0UL<rowindex>, 0UL<colindex>, 33) ])
+
+    let matrix = fromCoordinateList coo
+
+    let mul x y =
+        match x, y with
+        | None, None -> None
+        | None, Some(b) -> Some(b)
+        | Some(a), None -> Some(a)
+        | Some(a), Some(b) -> Some(a * b)
+
+    let result = reduceRows mul matrix
+    let vectorCoordinates = Vector.toCoordinateList result
+
+    let expected =
+        Vector.CoordinateList(1UL<Vector.dataLength>, [ (0UL<Vector.index>, 33) ])
+
+    Assert.Equal(expected, vectorCoordinates)
+
+
+[<Fact>]
+let ``let reduceCols sum on square power of two matrix`` () =
+    let coo =
+        CoordinateList(
+            2UL<nrows>,
+            2UL<ncols>,
+            [ (0UL<rowindex>, 0UL<colindex>, 12)
+              (0UL<rowindex>, 1UL<colindex>, 15)
+              (1UL<rowindex>, 0UL<colindex>, 17)
+              (1UL<rowindex>, 1UL<colindex>, 3) ]
+        )
+
+    let matrix = fromCoordinateList coo
+
+    let add x y =
+        match x, y with
+        | None, None -> None
+        | None, Some(b) -> Some(b)
+        | Some(a), None -> Some(a)
+        | Some(a), Some(b) -> Some(a + b)
+
+    let result = reduceCols add matrix
+    let vectorCoordinates = Vector.toCoordinateList result
+
+    let expected =
+        Vector.CoordinateList(2UL<Vector.dataLength>, [ (0UL<Vector.index>, 29); (1UL<Vector.index>, 18) ])
+
+    Assert.Equal(expected, vectorCoordinates)
+
+[<Fact>]
+let ``let reduceCols sum on square power of two matrix with empty col`` () =
+    let coo =
+        CoordinateList(
+            2UL<nrows>,
+            2UL<ncols>,
+            [ (0UL<rowindex>, 0UL<colindex>, 17); (1UL<rowindex>, 0UL<colindex>, 3) ]
+        )
+
+    let matrix = fromCoordinateList coo
+
+    let add x y =
+        match x, y with
+        | None, None -> None
+        | None, Some(b) -> Some(b)
+        | Some(a), None -> Some(a)
+        | Some(a), Some(b) -> Some(a + b)
+
+    let result = reduceCols add matrix
+    let vectorCoordinates = Vector.toCoordinateList result
+
+    let expected =
+        Vector.CoordinateList(2UL<Vector.dataLength>, [ (0UL<Vector.index>, 20) ])
+
+    Assert.Equal(expected, vectorCoordinates)
+
+[<Fact>]
+let ``let reduceCols sum on square not power of two matrix`` () =
+    let coo =
+        CoordinateList(
+            3UL<nrows>,
+            3UL<ncols>,
+            [ (0UL<rowindex>, 0UL<colindex>, 5)
+              (0UL<rowindex>, 1UL<colindex>, 7)
+              (0UL<rowindex>, 2UL<colindex>, 9)
+              (1UL<rowindex>, 1UL<colindex>, 11)
+              (1UL<rowindex>, 2UL<colindex>, 13)
+              (2UL<rowindex>, 0UL<colindex>, 15)
+              (2UL<rowindex>, 1UL<colindex>, 17) ]
+        )
+
+    let matrix = fromCoordinateList coo
+
+    let add x y =
+        match x, y with
+        | None, None -> None
+        | None, Some(b) -> Some(b)
+        | Some(a), None -> Some(a)
+        | Some(a), Some(b) -> Some(a + b)
+
+    let result = reduceCols add matrix
+    let vectorCoordinates = Vector.toCoordinateList result
+
+    let expected =
+        Vector.CoordinateList(
+            3UL<Vector.dataLength>,
+            [ (0UL<Vector.index>, 20); (1UL<Vector.index>, 35); (2UL<Vector.index>, 22) ]
+        )
+
+    Assert.Equal(expected, vectorCoordinates)
+
+[<Fact>]
+let ``let reduceCols mul on square not power of two matrix`` () =
+    let coo =
+        CoordinateList(
+            3UL<nrows>,
+            3UL<ncols>,
+            [ (0UL<rowindex>, 0UL<colindex>, 5)
+              (0UL<rowindex>, 1UL<colindex>, 7)
+              (0UL<rowindex>, 2UL<colindex>, 9)
+              (1UL<rowindex>, 1UL<colindex>, 11)
+              (1UL<rowindex>, 2UL<colindex>, 13)
+              (2UL<rowindex>, 0UL<colindex>, 15)
+              (2UL<rowindex>, 1UL<colindex>, 17) ]
+        )
+
+    let matrix = fromCoordinateList coo
+
+    let mul x y =
+        match x, y with
+        | None, None -> None
+        | None, Some(b) -> Some(b)
+        | Some(a), None -> Some(a)
+        | Some(a), Some(b) -> Some(a * b)
+
+    let result = reduceCols mul matrix
+    let vectorCoordinates = Vector.toCoordinateList result
+
+    let expected =
+        Vector.CoordinateList(
+            3UL<Vector.dataLength>,
+            [ (0UL<Vector.index>, 75); (1UL<Vector.index>, 1309); (2UL<Vector.index>, 117) ]
+        )
+
+    Assert.Equal(expected, vectorCoordinates)
+
+[<Fact>]
+let ``let reduceCols sum on rectangular matrix`` () =
+    let coo =
+        CoordinateList(
+            2UL<nrows>,
+            3UL<ncols>,
+            [ (0UL<rowindex>, 0UL<colindex>, 5)
+              (0UL<rowindex>, 1UL<colindex>, 7)
+              (0UL<rowindex>, 2UL<colindex>, 9)
+              (1UL<rowindex>, 1UL<colindex>, 11)
+              (1UL<rowindex>, 2UL<colindex>, 13) ]
+        )
+
+    let matrix = fromCoordinateList coo
+
+    let sum x y =
+        match x, y with
+        | None, None -> None
+        | None, Some(b) -> Some(b)
+        | Some(a), None -> Some(a)
+        | Some(a), Some(b) -> Some(a + b)
+
+    let result = reduceCols sum matrix
+    let vectorCoordinates = Vector.toCoordinateList result
+
+    let expected =
+        Vector.CoordinateList(
+            3UL<Vector.dataLength>,
+            [ (0UL<Vector.index>, 5); (1UL<Vector.index>, 18); (2UL<Vector.index>, 22) ]
+        )
+
+    Assert.Equal(expected, vectorCoordinates)
+
+[<Fact>]
+let ``let reduceCols sum on empty matrix`` () =
+    let coo = CoordinateList(2UL<nrows>, 3UL<ncols>, [])
+    let matrix = fromCoordinateList coo
+
+    let sum x y =
+        match x, y with
+        | None, None -> None
+        | None, Some(b) -> Some(b)
+        | Some(a), None -> Some(a)
+        | Some(a), Some(b) -> Some(a + b)
+
+    let result = reduceCols sum matrix
+    let vectorCoordinates = Vector.toCoordinateList result
+    let expected = Vector.CoordinateList(3UL<Vector.dataLength>, [])
+    Assert.Equal(expected, vectorCoordinates)
+
+[<Fact>]
+let ``let reduceCols mul on single matrix`` () =
+    let coo =
+        CoordinateList(1UL<nrows>, 1UL<ncols>, [ (0UL<rowindex>, 0UL<colindex>, 33) ])
+
+    let matrix = fromCoordinateList coo
+
+    let mul x y =
+        match x, y with
+        | None, None -> None
+        | None, Some(b) -> Some(b)
+        | Some(a), None -> Some(a)
+        | Some(a), Some(b) -> Some(a * b)
+
+    let result = reduceCols mul matrix
+    let vectorCoordinates = Vector.toCoordinateList result
+
+    let expected =
+        Vector.CoordinateList(1UL<Vector.dataLength>, [ (0UL<Vector.index>, 33) ])
+
+    Assert.Equal(expected, vectorCoordinates)
+
+[<Fact>]
+let ``kronecker product with square power of 2 x square power of two matrixes`` () =
+    let cooA =
+        CoordinateList(
+            2UL<nrows>,
+            2UL<ncols>,
+            [ (0UL<rowindex>, 0UL<colindex>, 1)
+              (0UL<rowindex>, 1UL<colindex>, 2)
+              (1UL<rowindex>, 0UL<colindex>, 3)
+              (1UL<rowindex>, 1UL<colindex>, 4) ]
+        )
+
+    let A = fromCoordinateList cooA
+
+    let cooB =
+        CoordinateList(
+            2UL<nrows>,
+            2UL<ncols>,
+            [ (0UL<rowindex>, 0UL<colindex>, 5)
+              (0UL<rowindex>, 1UL<colindex>, 6)
+              (1UL<rowindex>, 0UL<colindex>, 7)
+              (1UL<rowindex>, 1UL<colindex>, 8) ]
+        )
+
+    let B = fromCoordinateList cooB
+
+    let multiplyOp a b = Some(a * b)
+
+    match kroneckerProduct A B multiplyOp with
+    | Error msg -> Assert.True(false, msg)
+    | Ok result ->
+        let coo = toCoordinateList result
+
+        let expected =
+            CoordinateList(
+                4UL<nrows>,
+                4UL<ncols>,
+                [ (0UL<rowindex>, 0UL<colindex>, 5)
+                  (0UL<rowindex>, 1UL<colindex>, 6)
+                  (0UL<rowindex>, 2UL<colindex>, 10)
+                  (0UL<rowindex>, 3UL<colindex>, 12)
+                  (1UL<rowindex>, 0UL<colindex>, 7)
+                  (1UL<rowindex>, 1UL<colindex>, 8)
+                  (1UL<rowindex>, 2UL<colindex>, 14)
+                  (1UL<rowindex>, 3UL<colindex>, 16)
+                  (2UL<rowindex>, 0UL<colindex>, 15)
+                  (2UL<rowindex>, 1UL<colindex>, 18)
+                  (2UL<rowindex>, 2UL<colindex>, 20)
+                  (2UL<rowindex>, 3UL<colindex>, 24)
+                  (3UL<rowindex>, 0UL<colindex>, 21)
+                  (3UL<rowindex>, 1UL<colindex>, 24)
+                  (3UL<rowindex>, 2UL<colindex>, 28)
+                  (3UL<rowindex>, 3UL<colindex>, 32) ]
+            )
+
+        Assert.Equal(expected, coo)
+
+[<Fact>]
+let ``kronecker product with square not power of 2 x square not power of two matrixes`` () =
+    let cooA =
+        CoordinateList(
+            3UL<nrows>,
+            3UL<ncols>,
+            [ (0UL<rowindex>, 0UL<colindex>, 1)
+              (0UL<rowindex>, 1UL<colindex>, 2)
+              (0UL<rowindex>, 2UL<colindex>, 3)
+              (1UL<rowindex>, 0UL<colindex>, 4)
+              (1UL<rowindex>, 1UL<colindex>, 5)
+              (1UL<rowindex>, 2UL<colindex>, 6)
+              (2UL<rowindex>, 0UL<colindex>, 7)
+              (2UL<rowindex>, 1UL<colindex>, 8)
+              (2UL<rowindex>, 2UL<colindex>, 9) ]
+        )
+
+    let A = fromCoordinateList cooA
+
+    let cooB =
+        CoordinateList(
+            3UL<nrows>,
+            3UL<ncols>,
+            [ (0UL<rowindex>, 0UL<colindex>, 10)
+              (0UL<rowindex>, 1UL<colindex>, 11)
+              (0UL<rowindex>, 2UL<colindex>, 12)
+              (1UL<rowindex>, 0UL<colindex>, 13)
+              (1UL<rowindex>, 1UL<colindex>, 14)
+              (1UL<rowindex>, 2UL<colindex>, 15)
+              (2UL<rowindex>, 0UL<colindex>, 16)
+              (2UL<rowindex>, 1UL<colindex>, 17)
+              (2UL<rowindex>, 2UL<colindex>, 18) ]
+        )
+
+    let B = fromCoordinateList cooB
+
+    let multiplyOp a b = Some(a * b)
+
+    match kroneckerProduct A B multiplyOp with
+    | Error msg -> Assert.True(false, msg)
+    | Ok result ->
+        let coo = toCoordinateList result
+
+        let expected =
+            CoordinateList(
+                9UL<nrows>,
+                9UL<ncols>,
+                [ (0UL<rowindex>, 0UL<colindex>, 10)
+                  (0UL<rowindex>, 1UL<colindex>, 11)
+                  (0UL<rowindex>, 2UL<colindex>, 12)
+                  (0UL<rowindex>, 3UL<colindex>, 20)
+                  (0UL<rowindex>, 4UL<colindex>, 22)
+                  (0UL<rowindex>, 5UL<colindex>, 24)
+                  (0UL<rowindex>, 6UL<colindex>, 30)
+                  (0UL<rowindex>, 7UL<colindex>, 33)
+                  (0UL<rowindex>, 8UL<colindex>, 36)
+                  (1UL<rowindex>, 0UL<colindex>, 13)
+                  (1UL<rowindex>, 1UL<colindex>, 14)
+                  (1UL<rowindex>, 2UL<colindex>, 15)
+                  (1UL<rowindex>, 3UL<colindex>, 26)
+                  (1UL<rowindex>, 4UL<colindex>, 28)
+                  (1UL<rowindex>, 5UL<colindex>, 30)
+                  (1UL<rowindex>, 6UL<colindex>, 39)
+                  (1UL<rowindex>, 7UL<colindex>, 42)
+                  (1UL<rowindex>, 8UL<colindex>, 45)
+                  (2UL<rowindex>, 0UL<colindex>, 16)
+                  (2UL<rowindex>, 1UL<colindex>, 17)
+                  (2UL<rowindex>, 2UL<colindex>, 18)
+                  (2UL<rowindex>, 3UL<colindex>, 32)
+                  (2UL<rowindex>, 4UL<colindex>, 34)
+                  (2UL<rowindex>, 5UL<colindex>, 36)
+                  (2UL<rowindex>, 6UL<colindex>, 48)
+                  (2UL<rowindex>, 7UL<colindex>, 51)
+                  (2UL<rowindex>, 8UL<colindex>, 54)
+                  (3UL<rowindex>, 0UL<colindex>, 40)
+                  (3UL<rowindex>, 1UL<colindex>, 44)
+                  (3UL<rowindex>, 2UL<colindex>, 48)
+                  (3UL<rowindex>, 3UL<colindex>, 50)
+                  (3UL<rowindex>, 4UL<colindex>, 55)
+                  (3UL<rowindex>, 5UL<colindex>, 60)
+                  (3UL<rowindex>, 6UL<colindex>, 60)
+                  (3UL<rowindex>, 7UL<colindex>, 66)
+                  (3UL<rowindex>, 8UL<colindex>, 72)
+                  (4UL<rowindex>, 0UL<colindex>, 52)
+                  (4UL<rowindex>, 1UL<colindex>, 56)
+                  (4UL<rowindex>, 2UL<colindex>, 60)
+                  (4UL<rowindex>, 3UL<colindex>, 65)
+                  (4UL<rowindex>, 4UL<colindex>, 70)
+                  (4UL<rowindex>, 5UL<colindex>, 75)
+                  (4UL<rowindex>, 6UL<colindex>, 78)
+                  (4UL<rowindex>, 7UL<colindex>, 84)
+                  (4UL<rowindex>, 8UL<colindex>, 90)
+                  (5UL<rowindex>, 0UL<colindex>, 64)
+                  (5UL<rowindex>, 1UL<colindex>, 68)
+                  (5UL<rowindex>, 2UL<colindex>, 72)
+                  (5UL<rowindex>, 3UL<colindex>, 80)
+                  (5UL<rowindex>, 4UL<colindex>, 85)
+                  (5UL<rowindex>, 5UL<colindex>, 90)
+                  (5UL<rowindex>, 6UL<colindex>, 96)
+                  (5UL<rowindex>, 7UL<colindex>, 102)
+                  (5UL<rowindex>, 8UL<colindex>, 108)
+                  (6UL<rowindex>, 0UL<colindex>, 70)
+                  (6UL<rowindex>, 1UL<colindex>, 77)
+                  (6UL<rowindex>, 2UL<colindex>, 84)
+                  (6UL<rowindex>, 3UL<colindex>, 80)
+                  (6UL<rowindex>, 4UL<colindex>, 88)
+                  (6UL<rowindex>, 5UL<colindex>, 96)
+                  (6UL<rowindex>, 6UL<colindex>, 90)
+                  (6UL<rowindex>, 7UL<colindex>, 99)
+                  (6UL<rowindex>, 8UL<colindex>, 108)
+                  (7UL<rowindex>, 0UL<colindex>, 91)
+                  (7UL<rowindex>, 1UL<colindex>, 98)
+                  (7UL<rowindex>, 2UL<colindex>, 105)
+                  (7UL<rowindex>, 3UL<colindex>, 104)
+                  (7UL<rowindex>, 4UL<colindex>, 112)
+                  (7UL<rowindex>, 5UL<colindex>, 120)
+                  (7UL<rowindex>, 6UL<colindex>, 117)
+                  (7UL<rowindex>, 7UL<colindex>, 126)
+                  (7UL<rowindex>, 8UL<colindex>, 135)
+                  (8UL<rowindex>, 0UL<colindex>, 112)
+                  (8UL<rowindex>, 1UL<colindex>, 119)
+                  (8UL<rowindex>, 2UL<colindex>, 126)
+                  (8UL<rowindex>, 3UL<colindex>, 128)
+                  (8UL<rowindex>, 4UL<colindex>, 136)
+                  (8UL<rowindex>, 5UL<colindex>, 144)
+                  (8UL<rowindex>, 6UL<colindex>, 144)
+                  (8UL<rowindex>, 7UL<colindex>, 153)
+                  (8UL<rowindex>, 8UL<colindex>, 162) ]
+            )
+
+        Assert.Equal(expected, coo)
+
+[<Fact>]
+let ``kronecker product with rectangular and square matrixes`` () =
+    let cooA =
+        CoordinateList(
+            2UL<nrows>,
+            3UL<ncols>,
+            [ (0UL<rowindex>, 0UL<colindex>, 1)
+              (0UL<rowindex>, 1UL<colindex>, 2)
+              (1UL<rowindex>, 0UL<colindex>, 3)
+              (1UL<rowindex>, 2UL<colindex>, 4) ]
+        )
+
+    let A = fromCoordinateList cooA
+
+    let cooB =
+        CoordinateList(
+            2UL<nrows>,
+            2UL<ncols>,
+            [ (0UL<rowindex>, 0UL<colindex>, 5)
+              (0UL<rowindex>, 1UL<colindex>, 6)
+              (1UL<rowindex>, 0UL<colindex>, 7)
+              (1UL<rowindex>, 1UL<colindex>, 8) ]
+        )
+
+    let B = fromCoordinateList cooB
+
+    let result = kroneckerProduct A B (fun a b -> Some(a * b))
+
+    match result with
+    | Error msg -> Assert.True(false, msg)
+    | Ok res ->
+        let coo = toCoordinateList res
+
+        let expectedElements =
+            [ (0UL<rowindex>, 0UL<colindex>, 5)
+              (0UL<rowindex>, 1UL<colindex>, 6)
+              (1UL<rowindex>, 0UL<colindex>, 7)
+              (1UL<rowindex>, 1UL<colindex>, 8)
+
+              (0UL<rowindex>, 2UL<colindex>, 10)
+              (0UL<rowindex>, 3UL<colindex>, 12)
+              (1UL<rowindex>, 2UL<colindex>, 14)
+              (1UL<rowindex>, 3UL<colindex>, 16)
+
+              (2UL<rowindex>, 0UL<colindex>, 15)
+              (2UL<rowindex>, 1UL<colindex>, 18)
+              (3UL<rowindex>, 0UL<colindex>, 21)
+              (3UL<rowindex>, 1UL<colindex>, 24)
+
+              (2UL<rowindex>, 4UL<colindex>, 20)
+              (2UL<rowindex>, 5UL<colindex>, 24)
+              (3UL<rowindex>, 4UL<colindex>, 28)
+              (3UL<rowindex>, 5UL<colindex>, 32) ]
+            |> List.sortBy (fun (r, c, _) -> (r, c))
+
+        let actualElements = coo.list |> List.sortBy (fun (r, c, _) -> (r, c))
+
+        Assert.Equal(4UL<nrows>, coo.nrows)
+        Assert.Equal(6UL<ncols>, coo.ncols)
+        Assert.Equal<seq<uint64<rowindex> * uint64<colindex> * int>>(expectedElements, actualElements)
+
+[<Fact>]
+let ``kronecker product with square and rectangular matrixes`` () =
+    let cooA =
+        CoordinateList(
+            2UL<nrows>,
+            2UL<ncols>,
+            [ (0UL<rowindex>, 0UL<colindex>, 1)
+              (0UL<rowindex>, 1UL<colindex>, 2)
+              (1UL<rowindex>, 0UL<colindex>, 3)
+              (1UL<rowindex>, 1UL<colindex>, 4) ]
+        )
+
+    let A = fromCoordinateList cooA
+
+    let cooB =
+        CoordinateList(
+            2UL<nrows>,
+            3UL<ncols>,
+            [ (0UL<rowindex>, 0UL<colindex>, 5)
+              (0UL<rowindex>, 1UL<colindex>, 6)
+              (0UL<rowindex>, 2UL<colindex>, 7)
+              (1UL<rowindex>, 0UL<colindex>, 8)
+              (1UL<rowindex>, 1UL<colindex>, 9)
+              (1UL<rowindex>, 2UL<colindex>, 10) ]
+        )
+
+    let B = fromCoordinateList cooB
+
+    let result = kroneckerProduct A B (fun a b -> Some(a * b))
+
+    match result with
+    | Error msg -> Assert.True(false, msg)
+    | Ok res ->
+        let coo = toCoordinateList res
+
+        let expectedElements =
+            [ (0UL<rowindex>, 0UL<colindex>, 5)
+              (0UL<rowindex>, 1UL<colindex>, 6)
+              (0UL<rowindex>, 2UL<colindex>, 7)
+              (1UL<rowindex>, 0UL<colindex>, 8)
+              (1UL<rowindex>, 1UL<colindex>, 9)
+              (1UL<rowindex>, 2UL<colindex>, 10)
+
+              (0UL<rowindex>, 3UL<colindex>, 10)
+              (0UL<rowindex>, 4UL<colindex>, 12)
+              (0UL<rowindex>, 5UL<colindex>, 14)
+              (1UL<rowindex>, 3UL<colindex>, 16)
+              (1UL<rowindex>, 4UL<colindex>, 18)
+              (1UL<rowindex>, 5UL<colindex>, 20)
+
+              (2UL<rowindex>, 0UL<colindex>, 15)
+              (2UL<rowindex>, 1UL<colindex>, 18)
+              (2UL<rowindex>, 2UL<colindex>, 21)
+              (3UL<rowindex>, 0UL<colindex>, 24)
+              (3UL<rowindex>, 1UL<colindex>, 27)
+              (3UL<rowindex>, 2UL<colindex>, 30)
+
+              (2UL<rowindex>, 3UL<colindex>, 20)
+              (2UL<rowindex>, 4UL<colindex>, 24)
+              (2UL<rowindex>, 5UL<colindex>, 28)
+              (3UL<rowindex>, 3UL<colindex>, 32)
+              (3UL<rowindex>, 4UL<colindex>, 36)
+              (3UL<rowindex>, 5UL<colindex>, 40) ]
+            |> List.sortBy (fun (r, c, _) -> (r, c))
+
+        let actualElements = coo.list |> List.sortBy (fun (r, c, _) -> (r, c))
+
+        Assert.Equal(4UL<nrows>, coo.nrows)
+        Assert.Equal(6UL<ncols>, coo.ncols)
+        Assert.Equal<seq<uint64<rowindex> * uint64<colindex> * int>>(expectedElements, actualElements)
+
+[<Fact>]
+let ``kronecker product of matrix with empty matrix`` () =
+    let cooA =
+        CoordinateList(
+            2UL<nrows>,
+            2UL<ncols>,
+            [ (0UL<rowindex>, 0UL<colindex>, 1)
+              (0UL<rowindex>, 1UL<colindex>, 2)
+              (1UL<rowindex>, 0UL<colindex>, 3)
+              (1UL<rowindex>, 1UL<colindex>, 4) ]
+        )
+
+    let A = fromCoordinateList cooA
+
+    let emptyMatrix = fromCoordinateList (CoordinateList(2UL<nrows>, 2UL<ncols>, []))
+
+    let result = kroneckerProduct A emptyMatrix (fun a b -> Some(a * b))
+
+    match result with
+    | Error msg -> Assert.True(false, msg)
+    | Ok res ->
+        let coo = toCoordinateList res
+        let expected = CoordinateList(4UL<nrows>, 4UL<ncols>, [])
+        Assert.Equal(expected, coo)
+
+[<Fact>]
+let ``kronecker product of empty matrix with matrix`` () =
+    let emptyMatrix = fromCoordinateList (CoordinateList(2UL<nrows>, 2UL<ncols>, []))
+
+    let cooB =
+        CoordinateList(
+            2UL<nrows>,
+            2UL<ncols>,
+            [ (0UL<rowindex>, 0UL<colindex>, 1)
+              (0UL<rowindex>, 1UL<colindex>, 2)
+              (1UL<rowindex>, 0UL<colindex>, 3)
+              (1UL<rowindex>, 1UL<colindex>, 4) ]
+        )
+
+    let B = fromCoordinateList cooB
+
+    let result = kroneckerProduct emptyMatrix B (fun a b -> Some(a * b))
+
+    match result with
+    | Error msg -> Assert.True(false, msg)
+    | Ok res ->
+        let coo = toCoordinateList res
+        let expected = CoordinateList(4UL<nrows>, 4UL<ncols>, [])
+        Assert.Equal(expected, coo)
+
+[<Fact>]
+let ``kronecker product of matrix with zeros`` () =
+    let cooA =
+        CoordinateList(1UL<nrows>, 1UL<ncols>, [ (0UL<rowindex>, 0UL<colindex>, 2) ])
+
+    let A = fromCoordinateList cooA
+
+    let cooB =
+        CoordinateList(2UL<nrows>, 2UL<ncols>, [ (0UL<rowindex>, 0UL<colindex>, 0); (1UL<rowindex>, 1UL<colindex>, 3) ])
+
+    let B = fromCoordinateList cooB
+
+    let result = kroneckerProduct A B (fun a b -> Some(a * b))
+
+    match result with
+    | Error msg -> Assert.True(false, msg)
+    | Ok res ->
+        let coo = toCoordinateList res
+
+        let expectedElements =
+            [ (0UL<rowindex>, 0UL<colindex>, 0); (1UL<rowindex>, 1UL<colindex>, 6) ]
+            |> List.sortBy (fun (r, c, _) -> (r, c))
+
+        let actualElements = coo.list |> List.sortBy (fun (r, c, _) -> (r, c))
+
+        Assert.Equal(2UL<nrows>, coo.nrows)
+        Assert.Equal(2UL<ncols>, coo.ncols)
+        Assert.Equal<seq<uint64<rowindex> * uint64<colindex> * int>>(expectedElements, actualElements)
+
+[<Fact>]
+let ``kronecker product resulting entirely in explicit zeros`` () =
+    let cooA =
+        CoordinateList(2UL<nrows>, 2UL<ncols>, [ (0UL<rowindex>, 0UL<colindex>, 0); (1UL<rowindex>, 1UL<colindex>, 0) ])
+
+    let A = fromCoordinateList cooA
+
+    let cooB =
+        CoordinateList(1UL<nrows>, 1UL<ncols>, [ (0UL<rowindex>, 0UL<colindex>, 5) ])
+
+    let B = fromCoordinateList cooB
+
+    let result = kroneckerProduct A B (fun a b -> Some(a * b))
+
+    match result with
+    | Error msg -> Assert.True(false, msg)
+    | Ok res ->
+        let coo = toCoordinateList res
+
+        let expectedElements =
+            [ (0UL<rowindex>, 0UL<colindex>, 0); (1UL<rowindex>, 1UL<colindex>, 0) ]
+            |> List.sortBy (fun (r, c, _) -> (r, c))
+
+        let actualElements = coo.list |> List.sortBy (fun (r, c, _) -> (r, c))
+
+        Assert.Equal(2UL<nrows>, coo.nrows)
+        Assert.Equal(2UL<ncols>, coo.ncols)
+        Assert.Equal<seq<uint64<rowindex> * uint64<colindex> * int>>(expectedElements, actualElements)
+
+[<Fact>]
+let ``kronecker product of square matrix with matrix 1x1`` () =
+    let cooA =
+        CoordinateList(
+            2UL<nrows>,
+            2UL<ncols>,
+            [ (0UL<rowindex>, 0UL<colindex>, 1)
+              (0UL<rowindex>, 1UL<colindex>, 2)
+              (1UL<rowindex>, 0UL<colindex>, 3)
+              (1UL<rowindex>, 1UL<colindex>, 4) ]
+        )
+
+    let A = fromCoordinateList cooA
+
+    let cooB =
+        CoordinateList(1UL<nrows>, 1UL<ncols>, [ (0UL<rowindex>, 0UL<colindex>, 3) ])
+
+    let B = fromCoordinateList cooB
+
+    let result = kroneckerProduct A B (fun a b -> Some(a * b))
+
+    match result with
+    | Error msg -> Assert.True(false, msg)
+    | Ok res ->
+        let coo = toCoordinateList res
+
+        let expected =
+            CoordinateList(
+                2UL<nrows>,
+                2UL<ncols>,
+                [ (0UL<rowindex>, 0UL<colindex>, 3)
+                  (0UL<rowindex>, 1UL<colindex>, 6)
+                  (1UL<rowindex>, 0UL<colindex>, 9)
+                  (1UL<rowindex>, 1UL<colindex>, 12) ]
+            )
+
+        Assert.Equal(expected, coo)
+
+[<Fact>]
+let ``kronecker product of matrix 1x1 with square matrix`` () =
+    let cooA =
+        CoordinateList(1UL<nrows>, 1UL<ncols>, [ (0UL<rowindex>, 0UL<colindex>, 3) ])
+
+    let A = fromCoordinateList cooA
+
+    let cooB =
+        CoordinateList(
+            2UL<nrows>,
+            2UL<ncols>,
+            [ (0UL<rowindex>, 0UL<colindex>, 1)
+              (0UL<rowindex>, 1UL<colindex>, 2)
+              (1UL<rowindex>, 0UL<colindex>, 3)
+              (1UL<rowindex>, 1UL<colindex>, 4) ]
+        )
+
+    let B = fromCoordinateList cooB
+
+    let result = kroneckerProduct A B (fun a b -> Some(a * b))
+
+    match result with
+    | Error msg -> Assert.True(false, msg)
+    | Ok res ->
+        let coo = toCoordinateList res
+
+        let expected =
+            CoordinateList(
+                2UL<nrows>,
+                2UL<ncols>,
+                [ (0UL<rowindex>, 0UL<colindex>, 3)
+                  (0UL<rowindex>, 1UL<colindex>, 6)
+                  (1UL<rowindex>, 0UL<colindex>, 9)
+                  (1UL<rowindex>, 1UL<colindex>, 12) ]
+            )
+
+        Assert.Equal(expected, coo)
+
+[<Fact>]
+let ``kronecker dimension check`` () =
+    let cooA =
+        CoordinateList(3UL<nrows>, 4UL<ncols>, [ (0UL<rowindex>, 0UL<colindex>, 1) ])
+
+    let A = fromCoordinateList cooA
+
+    let cooB =
+        CoordinateList(2UL<nrows>, 5UL<ncols>, [ (0UL<rowindex>, 0UL<colindex>, 1) ])
+
+    let B = fromCoordinateList cooB
+
+    let result = kroneckerProduct A B (fun a b -> Some(a * b))
+
+    match result with
+    | Error msg -> Assert.True(false, msg)
+    | Ok res ->
+        Assert.Equal(6UL<nrows>, res.nrows)
+        Assert.Equal(20UL<ncols>, res.ncols)
+
+[<Fact>]
+let ``kronecker product with sparse matrix on dense matrix`` () =
+    let cooA =
+        CoordinateList(10UL<nrows>, 10UL<ncols>, [ (5UL<rowindex>, 5UL<colindex>, 2) ])
+
+    let A = fromCoordinateList cooA
+
+    let cooB =
+        CoordinateList(
+            3UL<nrows>,
+            3UL<ncols>,
+            [ (0UL<rowindex>, 0UL<colindex>, 1)
+              (0UL<rowindex>, 1UL<colindex>, 2)
+              (0UL<rowindex>, 2UL<colindex>, 3)
+              (1UL<rowindex>, 0UL<colindex>, 4)
+              (1UL<rowindex>, 1UL<colindex>, 5)
+              (1UL<rowindex>, 2UL<colindex>, 6)
+              (2UL<rowindex>, 0UL<colindex>, 7)
+              (2UL<rowindex>, 1UL<colindex>, 8)
+              (2UL<rowindex>, 2UL<colindex>, 9) ]
+        )
+
+    let B = fromCoordinateList cooB
+
+    let result = kroneckerProduct A B (fun a b -> Some(a * b))
+
+    match result with
+    | Error msg -> Assert.True(false, msg)
+    | Ok res ->
+        let coo = toCoordinateList res
+
+        let expected =
+            CoordinateList(
+                30UL<nrows>,
+                30UL<ncols>,
+                [ (15UL<rowindex>, 15UL<colindex>, 2)
+                  (15UL<rowindex>, 16UL<colindex>, 4)
+                  (15UL<rowindex>, 17UL<colindex>, 6)
+                  (16UL<rowindex>, 15UL<colindex>, 8)
+                  (16UL<rowindex>, 16UL<colindex>, 10)
+                  (16UL<rowindex>, 17UL<colindex>, 12)
+                  (17UL<rowindex>, 15UL<colindex>, 14)
+                  (17UL<rowindex>, 16UL<colindex>, 16)
+                  (17UL<rowindex>, 17UL<colindex>, 18) ]
+            )
+
+        Assert.Equal(expected, coo)
+
+[<Fact>]
+let ``kronecker product with filtering (only even results)`` () =
+    let cooA =
+        CoordinateList(
+            2UL<nrows>,
+            2UL<ncols>,
+            [ (0UL<rowindex>, 0UL<colindex>, 1)
+              (0UL<rowindex>, 1UL<colindex>, 2)
+              (1UL<rowindex>, 0UL<colindex>, 3)
+              (1UL<rowindex>, 1UL<colindex>, 4) ]
+        )
+
+    let A = fromCoordinateList cooA
+
+    let cooB =
+        CoordinateList(
+            2UL<nrows>,
+            2UL<ncols>,
+            [ (0UL<rowindex>, 0UL<colindex>, 1)
+              (0UL<rowindex>, 1UL<colindex>, 2)
+              (1UL<rowindex>, 0UL<colindex>, 3)
+              (1UL<rowindex>, 1UL<colindex>, 4) ]
+        )
+
+    let B = fromCoordinateList cooB
+
+    let evenOnly a b =
+        let prod = a * b
+        if prod % 2 = 0 then Some prod else None
+
+    let result = kroneckerProduct A B evenOnly
+
+    match result with
+    | Error msg -> Assert.True(false, msg)
+    | Ok res ->
+        let coo = toCoordinateList res
+        Assert.True(coo.list |> List.forall (fun (_, _, v) -> v % 2 = 0))
+        Assert.True(List.length coo.list < 16)
diff --git a/QuadTree.Tests/Tests.Vector.fs b/QuadTree.Tests/Tests.Vector.fs
index 4eb1af9..61ee41f 100644
--- a/QuadTree.Tests/Tests.Vector.fs
+++ b/QuadTree.Tests/Tests.Vector.fs
@@ -1,6 +1,7 @@
 module Vector.Tests
 
 open Xunit
+open System
 
 open Vector
 open Common
@@ -903,3 +904,252 @@ let ``Init vector`` () =
 
     let actual = Vector.init 3UL<dataLength> (fun i -> Some(int i))
     Assert.Equal(expected, actual)
+
+
+[<Fact>]
+let ``map on empty vector returns empty vector`` () =
+    let vec = fromCoordinateList (CoordinateList(0UL<dataLength>, []))
+
+    let result =
+        map vec (fun (x: int option) ->
+            match x with
+            | Some v -> Some(v * 3)
+            | None -> None)
+
+    let coo = toCoordinateList result
+    Assert.Equal(CoordinateList(0UL<dataLength>, []), coo)
+
+[<Fact>]
+let ``map with function that turns all the elements into zeros`` () =
+    let vec =
+        fromCoordinateList (CoordinateList(7UL<dataLength>, [ (1UL<index>, 3); (2UL<index>, 3); (5UL<index>, 100) ]))
+
+    let result = map vec (fun _ -> Some 0)
+    let coo = toCoordinateList result
+
+    Assert.Equal(
+        CoordinateList(
+            7UL<dataLength>,
+            [ (0UL<index>, 0)
+              (1UL<index>, 0)
+              (2UL<index>, 0)
+              (3UL<index>, 0)
+              (4UL<index>, 0)
+              (5UL<index>, 0)
+              (6UL<index>, 0) ]
+        ),
+        coo
+    )
+
+[<Fact>]
+let ``map with function that resets all the elements`` () =
+    let vec =
+        fromCoordinateList (CoordinateList(7UL<dataLength>, [ (1UL<index>, 3); (2UL<index>, 3); (5UL<index>, 100) ]))
+
+    let result = map vec (fun _ -> None)
+    let coo = toCoordinateList result
+    Assert.Equal(CoordinateList(7UL<dataLength>, []), coo)
+
+[<Fact>]
+let ``map can change type from int to string`` () =
+    let vec =
+        fromCoordinateList (CoordinateList(3UL<dataLength>, [ (1UL<index>, 11); (2UL<index>, 33) ]))
+
+    let result =
+        map vec (fun (x: int option) ->
+            match x with
+            | Some v -> Some(sprintf "str %d" v)
+            | None -> None)
+
+    let coo = toCoordinateList result
+    Assert.Equal(CoordinateList(3UL<dataLength>, [ (1UL<index>, "str 11"); (2UL<index>, "str 33") ]), coo)
+
+[<Fact>]
+let ``fromCoordinateList with out-of-range index throws exception`` () =
+    let coo = CoordinateList(6UL<dataLength>, [ (9UL<index>, 8) ])
+    Assert.Throws<Exception>(fun () -> fromCoordinateList coo |> ignore)
+
+[<Fact>]
+let ``fromCoordinateList with unsorted coordinates works correctly`` () =
+    let coo =
+        CoordinateList(7UL<dataLength>, [ (5UL<index>, 3); (3UL<index>, 2); (1UL<index>, 100) ])
+
+    let vec = fromCoordinateList coo
+    let result = toCoordinateList vec
+    Assert.Equal(CoordinateList(7UL<dataLength>, [ (1UL<index>, 100); (3UL<index>, 2); (5UL<index>, 3) ]), result)
+
+[<Fact>]
+let ``fromCoordinateList with duplicate indicies returns the last of them`` () =
+    let coo = CoordinateList(3UL<dataLength>, [ (1UL<index>, 33); (1UL<index>, 100) ])
+    let vec = Vector.fromCoordinateList coo
+    let result = Vector.toCoordinateList vec
+    Assert.Equal(CoordinateList(3UL<dataLength>, [ (1UL<index>, 100) ]), result)
+
+[<Fact>]
+let ``fromCoordinateList with zero size throws Exception`` () =
+    let coo = CoordinateList(0UL<dataLength>, [ (33UL<index>, 33); (39UL<index>, 1) ])
+    Assert.Throws<Exception>(fun () -> fromCoordinateList coo |> ignore)
+
+[<Fact>]
+let ``mapi works with index`` () =
+    let vec =
+        fromCoordinateList (CoordinateList(5UL<dataLength>, [ (2UL<index>, 10); (3UL<index>, 23) ]))
+
+    let result =
+        mapi vec (fun (i: uint64<index>) (x: int option) ->
+            match x with
+            | Some v -> Some(int i + v)
+            | None -> None)
+
+    let coo = toCoordinateList result
+    Assert.Equal(CoordinateList(5UL<dataLength>, [ (2UL<index>, 12); (3UL<index>, 26) ]), coo)
+
+[<Fact>]
+let ``mapi on empty vector returns empty vector`` () =
+    let vec = fromCoordinateList (CoordinateList(0UL<dataLength>, []))
+
+    let result =
+        mapi vec (fun (i: uint64<index>) (x: int option) ->
+            match x with
+            | Some v -> Some(int i + v)
+            | None -> None)
+
+    let coo = toCoordinateList result
+    Assert.Equal(CoordinateList(0UL<dataLength>, []), coo)
+
+[<Fact>]
+let ``mapi with special function returns empty vector`` () =
+    let vec =
+        fromCoordinateList (CoordinateList(5UL<dataLength>, [ (1UL<index>, 33); (3UL<index>, 88) ]))
+
+    let result =
+        mapi vec (fun (i: uint64<index>) (x: int option) ->
+            match x with
+            | Some v ->
+                match int i % 2 with
+                | 0 -> Some(v * 2)
+                | _ -> None
+            | None -> None)
+
+    let coo = toCoordinateList result
+    Assert.Equal(CoordinateList(5UL<dataLength>, []), coo)
+
+[<Fact>]
+let ``mapi works with function does not depend on the index`` () =
+    let vec =
+        fromCoordinateList (CoordinateList(5UL<dataLength>, [ (2UL<index>, 10); (3UL<index>, 23) ]))
+
+    let result =
+        mapi vec (fun (i: uint64<index>) (x: int option) ->
+            match x with
+            | Some v -> Some(v + 4)
+            | None -> None)
+
+    let coo = toCoordinateList result
+    Assert.Equal(CoordinateList(5UL<dataLength>, [ (2UL<index>, 14); (3UL<index>, 27) ]), coo)
+
+[<Fact>]
+let ``slice returns error when start is negative`` () =
+    let vec = fromCoordinateList (CoordinateList(5UL<dataLength>, [ (2UL<index>, 10) ]))
+
+    match slice -1 3 vec with
+    | Result.Ok _ -> Assert.Fail("Expected Error")
+    | Result.Error msg -> Assert.Equal("Start should be >= 0", msg)
+
+[<Fact>]
+let ``slice returns error when end is negative`` () =
+    let vec = fromCoordinateList (CoordinateList(5UL<dataLength>, [ (2UL<index>, 10) ]))
+
+    match slice 1 -3 vec with
+    | Result.Ok _ -> Assert.Fail("Expected Error")
+    | Result.Error msg -> Assert.Equal("End should be >= 0", msg)
+
+[<Fact>]
+let ``slice returns error when start is out of range`` () =
+    let vec = fromCoordinateList (CoordinateList(5UL<dataLength>, [ (2UL<index>, 10) ]))
+
+    match slice 7 3 vec with
+    | Result.Ok _ -> Assert.Fail("Expected Error")
+    | Result.Error msg -> Assert.Equal("Start is out of Vector length", msg)
+
+[<Fact>]
+let ``slice returns error when end is out of range`` () =
+    let vec = fromCoordinateList (CoordinateList(5UL<dataLength>, [ (2UL<index>, 10) ]))
+
+    match slice 3 7 vec with
+    | Result.Ok _ -> Assert.Fail("Expected Error")
+    | Result.Error msg -> Assert.Equal("End is out of Vector length", msg)
+
+[<Fact>]
+let ``slice returns error when end is less than start`` () =
+    let vec = fromCoordinateList (CoordinateList(5UL<dataLength>, [ (2UL<index>, 10) ]))
+
+    match slice 4 3 vec with
+    | Result.Ok _ -> Assert.Fail("Expected Error")
+    | Result.Error msg -> Assert.Equal("End should be >= Start", msg)
+
+[<Fact>]
+let ``slice returns correct subvector`` () =
+    let vec =
+        fromCoordinateList (CoordinateList(7UL<dataLength>, [ (2UL<index>, 10); (6UL<index>, 20) ]))
+
+    match slice 1 3 vec with
+    | Result.Ok result ->
+        let coo = Vector.toCoordinateList result
+        Assert.Equal(CoordinateList(3UL<dataLength>, [ (1UL<index>, 10) ]), coo)
+    | Result.Error msg -> Assert.Fail(msg)
+
+[<Fact>]
+let ``slice returns subvector without elements`` () =
+    let vec =
+        fromCoordinateList (CoordinateList(7UL<dataLength>, [ (5UL<index>, 10); (6UL<index>, 20) ]))
+
+    match slice 0 2 vec with
+    | Result.Ok result ->
+        let coo = toCoordinateList result
+        Assert.Equal(CoordinateList(3UL<dataLength>, []), coo)
+    | Result.Error msg -> Assert.Fail(msg)
+
+[<Fact>]
+let ``slice returns single subvector`` () =
+    let vec =
+        fromCoordinateList (CoordinateList(7UL<dataLength>, [ (2UL<index>, 10); (6UL<index>, 20) ]))
+
+    match slice 2 2 vec with
+    | Result.Ok result ->
+        let coo = toCoordinateList result
+        Assert.Equal(CoordinateList(1UL<dataLength>, [ (0UL<index>, 10) ]), coo)
+    | Result.Error msg -> Assert.Fail(msg)
+
+[<Fact>]
+let ``slice returns correct subvector equals to vector`` () =
+    let vec =
+        fromCoordinateList (CoordinateList(7UL<dataLength>, [ (2UL<index>, 10); (3UL<index>, 33); (6UL<index>, 20) ]))
+
+    match slice 0 6 vec with
+    | Result.Ok result ->
+        let coo = toCoordinateList result
+        Assert.Equal(CoordinateList(7UL<dataLength>, [ (2UL<index>, 10); (3UL<index>, 33); (6UL<index>, 20) ]), coo)
+    | Result.Error msg -> Assert.Fail(msg)
+
+[<Fact>]
+let ``slice returns correct subvector when start of subvector equals to start of vector`` () =
+    let vec =
+        fromCoordinateList (CoordinateList(7UL<dataLength>, [ (0UL<index>, 10); (3UL<index>, 33); (6UL<index>, 20) ]))
+
+    match slice 0 4 vec with
+    | Result.Ok result ->
+        let coo = Vector.toCoordinateList result
+        Assert.Equal(CoordinateList(5UL<dataLength>, [ (0UL<index>, 10); (3UL<index>, 33) ]), coo)
+    | Result.Error msg -> Assert.Fail(msg)
+
+[<Fact>]
+let ``slice returns correct subvector when end of subvector equals to end of vector`` () =
+    let vec =
+        fromCoordinateList (CoordinateList(7UL<dataLength>, [ (0UL<index>, 10); (3UL<index>, 33); (6UL<index>, 20) ]))
+
+    match slice 2 6 vec with
+    | Result.Ok result ->
+        let coo = Vector.toCoordinateList result
+        Assert.Equal(CoordinateList(5UL<dataLength>, [ (1UL<index>, 33); (4UL<index>, 20) ]), coo)
+    | Result.Error msg -> Assert.Fail(msg)
diff --git a/QuadTree/Matrix.fs b/QuadTree/Matrix.fs
index 68ea7d7..e42917d 100644
--- a/QuadTree/Matrix.fs
+++ b/QuadTree/Matrix.fs
@@ -74,7 +74,20 @@ let private getQuadrantCoords (pr, pc) halfSize =
     (pr + halfSize * 1UL<rowindex>, pc + halfSize * 1UL<colindex>) // SOUTH EAST
 
 let fromCoordinateList (coo: CoordinateList<'a>) =
-    let nvals = (uint64 <| List.length coo.list) * 1UL<nvals>
+    let unique =
+        coo.list
+        |> List.groupBy (fun (i, j, _) -> (i, j))
+        |> List.map (fun ((i, j), entries) ->
+            let value = entries |> List.map (fun (_, _, v) -> v) |> List.last
+            (i, j, value))
+
+    if
+        unique
+        |> List.exists (fun (i, j, _) -> uint64 i >= uint64 coo.nrows || uint64 j >= uint64 coo.ncols)
+    then
+        failwith "Coordinates out of range"
+
+    let nvals = (uint64 <| List.length unique) * 1UL<nvals>
     let nrows = coo.nrows
     let ncols = coo.ncols
 
@@ -108,7 +121,7 @@ let fromCoordinateList (coo: CoordinateList<'a>) =
                 (traverse swCoo swp halfSize)
                 (traverse seCoo sep halfSize)
 
-    let tree = traverse coo.list (0UL<rowindex>, 0UL<colindex>) storageSize
+    let tree = traverse unique (0UL<rowindex>, 0UL<colindex>) storageSize
 
     SparseMatrix(nrows, ncols, nvals, Storage(storageSize * 1UL<storageSize>, tree))
 
@@ -135,11 +148,35 @@ let toCoordinateList (matrix: SparseMatrix<'a>) =
     let coo =
         traverse matrix.storage.data (0UL<rowindex>, 0UL<colindex>) (uint64 matrix.storage.size)
 
-    CoordinateList(nrows, ncols, coo)
+    let sorted = List.sort coo
+    CoordinateList(nrows, ncols, sorted)
 
 let empty nrows ncols =
     fromCoordinateList (CoordinateList(nrows, ncols, []))
 
+let map (matrix: SparseMatrix<'a>) f =
+    let rec inner (size: uint64<storageSize>) (tree: qtree<Option<'a>>) =
+        match tree with
+        | Node(nw, ne, sw, se) ->
+            let nwTree, nwNvals = inner (size / 2UL) nw
+            let neTree, neNvals = inner (size / 2UL) ne
+            let swTree, swNvals = inner (size / 2UL) sw
+            let seTree, seNvals = inner (size / 2UL) se
+            (mkNode nwTree neTree swTree seTree), nwNvals + neNvals + swNvals + seNvals
+        | Leaf(Dummy) -> Leaf(Dummy), 0UL<nvals>
+        | Leaf(UserValue(v)) ->
+            let res = f v
+
+            let nnz =
+                match res with
+                | None -> 0UL<nvals>
+                | _ -> (uint64 size) * (uint64 size) * 1UL<nvals>
+
+            Leaf(UserValue(res)), nnz
+
+    let newTree, newNvals = inner matrix.storage.size matrix.storage.data
+    SparseMatrix(matrix.nrows, matrix.ncols, newNvals, Storage(matrix.storage.size, newTree))
+
 let map2 (matrix1: SparseMatrix<_>) (matrix2: SparseMatrix<_>) f =
     let rec inner (size: uint64<storageSize>) matrix1 matrix2 =
         let _do x1 x2 x3 x4 y1 y2 y3 y4 =
@@ -394,3 +431,340 @@ let transpose (matrix: SparseMatrix<_>) =
 
 let mask (m1: SparseMatrix<'a>) (m2: SparseMatrix<'b>) f =
     map2 m1 m2 (fun m1 m2 -> if f m2 then m1 else None)
+
+let slice
+    (matrix: SparseMatrix<'a>)
+    (rowStart: int)
+    (rowEnd: int)
+    (colStart: int)
+    (colEnd: int)
+    : Result<SparseMatrix<'a>, string> =
+    if rowStart < 0 then
+        Error "Start row should be >= 0"
+    elif rowEnd < 0 then
+        Error "End row should be >= 0"
+    elif colStart < 0 then
+        Error "Start column should be >= 0"
+    elif colEnd < 0 then
+        Error "End column should be >= 0"
+    elif rowStart > int matrix.nrows - 1 then
+        Error "Start row is out of matrix length"
+    elif rowEnd > int matrix.nrows - 1 then
+        Error "End row is out of matrix length"
+    elif colStart > int matrix.ncols - 1 then
+        Error "Start column is out of matrix length"
+    elif colEnd > int matrix.ncols - 1 then
+        Error "End column is out of matrix length"
+    elif rowStart > rowEnd then
+        Error "Start row should be <= end row"
+    elif colStart > colEnd then
+        Error "Start column should be <= end column"
+    else
+        let rowStartIdx = uint64 rowStart * 1UL<rowindex>
+        let rowEndIdx = uint64 rowEnd * 1UL<rowindex>
+        let colStartIdx = uint64 colStart * 1UL<colindex>
+        let colEndIdx = uint64 colEnd * 1UL<colindex>
+        let newRows = uint64 (rowEnd - rowStart + 1) * 1UL<nrows>
+        let newCols = uint64 (colEnd - colStart + 1) * 1UL<ncols>
+
+        let newSize =
+            getNearestUpperPowerOfTwo (max (uint64 newRows) (uint64 newCols))
+            * 1UL<storageSize>
+
+        let rec cut (size: uint64<storageSize>) (row: uint64<rowindex>) (col: uint64<colindex>) tree =
+            let sizeRow = (uint64 size) * 1UL<rowindex>
+            let sizeCol = (uint64 size) * 1UL<colindex>
+
+            match tree with
+            | Node(nw, ne, sw, se) ->
+                let half = size / 2UL
+                let halfRow = (uint64 half) * 1UL<rowindex>
+                let halfCol = (uint64 half) * 1UL<colindex>
+
+                mkNode
+                    (cut half row col nw)
+                    (cut half row (col + halfCol) ne)
+                    (cut half (row + halfRow) col sw)
+                    (cut half (row + halfRow) (col + halfCol) se)
+            | Leaf(Dummy) -> Leaf Dummy
+            | Leaf(UserValue(v)) when
+                row >= rowStartIdx
+                && row + sizeRow - 1UL<rowindex> <= rowEndIdx
+                && col >= colStartIdx
+                && col + sizeCol - 1UL<colindex> <= colEndIdx
+                ->
+                Leaf(UserValue(v))
+            | _ -> Leaf Dummy
+
+        let cutTree =
+            cut matrix.storage.size 0UL<rowindex> 0UL<colindex> matrix.storage.data
+
+        let rec empty (size: uint64<storageSize>) =
+            match size with
+            | 1UL<storageSize> -> Leaf Dummy
+            | _ ->
+                let half = size / 2UL
+                let subtree = empty half
+                Node(subtree, subtree, subtree, subtree)
+
+        let rec insert (size: uint64<storageSize>) (row: uint64<rowindex>) (col: uint64<colindex>) value tree =
+            match size with
+            | 1UL<storageSize> -> Leaf(UserValue(value))
+            | _ ->
+                let half = size / 2UL
+                let halfRow = (uint64 half) * 1UL<rowindex>
+                let halfCol = (uint64 half) * 1UL<colindex>
+
+                match tree with
+                | Node(nw, ne, sw, se) ->
+                    if row < halfRow then
+                        if col < halfCol then
+                            Node(insert half row col value nw, ne, sw, se)
+                        else
+                            Node(nw, insert half row (col - halfCol) value ne, sw, se)
+                    else if col < halfCol then
+                        Node(nw, ne, insert half (row - halfRow) col value sw, se)
+                    else
+                        Node(nw, ne, sw, insert half (row - halfRow) (col - halfCol) value se)
+                | _ ->
+                    let emptySub = empty half
+
+                    if row < halfRow then
+                        if col < halfCol then
+                            Node(insert half row col value emptySub, emptySub, emptySub, emptySub)
+                        else
+                            Node(emptySub, insert half row (col - halfCol) value emptySub, emptySub, emptySub)
+                    else if col < halfCol then
+                        Node(emptySub, emptySub, insert half (row - halfRow) col value emptySub, emptySub)
+                    else
+                        Node(emptySub, emptySub, emptySub, insert half (row - halfRow) (col - halfCol) value emptySub)
+
+        let rec rebuild (size: uint64<storageSize>) (row: uint64<rowindex>) (col: uint64<colindex>) tree acc =
+            match tree with
+            | Leaf(Dummy) -> acc
+            | Leaf(UserValue(v)) ->
+                let newRow = row - uint64 rowStart * 1UL<rowindex>
+                let newCol = col - uint64 colStart * 1UL<colindex>
+                insert newSize newRow newCol v acc
+            | Node(nw, ne, sw, se) ->
+                let half = size / 2UL
+                let halfRow = (uint64 half) * 1UL<rowindex>
+                let halfCol = (uint64 half) * 1UL<colindex>
+                let acc = rebuild half row col nw acc
+                let acc = rebuild half row (col + halfCol) ne acc
+                let acc = rebuild half (row + halfRow) col sw acc
+                rebuild half (row + halfRow) (col + halfCol) se acc
+
+        let emptyTree = empty newSize
+
+        let shiftedTree =
+            rebuild matrix.storage.size 0UL<rowindex> 0UL<colindex> cutTree emptyTree
+
+        let rec count (size: uint64<storageSize>) tree =
+            match tree with
+            | Node(nw, ne, sw, se) ->
+                let half = size / 2UL
+                count half nw + count half ne + count half sw + count half se
+            | Leaf(UserValue(_)) -> 1UL<nvals>
+            | _ -> 0UL<nvals>
+
+        let nvals = count newSize shiftedTree
+
+        Ok(SparseMatrix(newRows, newCols, nvals, Storage(newSize, shiftedTree)))
+
+let reduceRows (op: 'a option -> 'a option -> 'a option) (matrix: SparseMatrix<'a>) : Vector.SparseVector<'a> =
+    let rows = matrix.nrows
+
+    let rec inner (size: uint64<storageSize>) (row: uint64<rowindex>) (col: uint64<colindex>) tree acc =
+        match tree with
+        | Node(nw, ne, sw, se) ->
+            let half = size / 2UL
+            let halfRow = (uint64 half) * 1UL<rowindex>
+            let halfCol = (uint64 half) * 1UL<colindex>
+            let acc = inner half row col nw acc
+            let acc = inner half row (col + halfCol) ne acc
+            let acc = inner half (row + halfRow) col sw acc
+            let acc = inner half (row + halfRow) (col + halfCol) se acc
+            acc
+        | Leaf(Dummy) -> acc
+        | Leaf(UserValue(None)) -> acc
+        | Leaf(UserValue(Some v)) -> (row, v) :: acc
+
+    let pairs =
+        inner matrix.storage.size 0UL<rowindex> 0UL<colindex> matrix.storage.data []
+
+    let grouped =
+        List.sortBy fst pairs
+        |> List.groupBy fst
+        |> List.map (fun (row, rowSet) -> row, rowSet |> List.map snd |> List.map Some)
+
+    let reduced =
+        grouped
+        |> List.map (fun (row, values) ->
+            match values with
+            | [] -> row, None
+            | head :: tail -> row, List.fold op head tail)
+
+    let data =
+        reduced
+        |> List.choose (fun (row, v) ->
+            match v with
+            | None -> None
+            | Some x -> Some(row, 0UL<colindex>, x))
+        |> List.sort
+
+    let vectorData =
+        data |> List.map (fun (row, _, v) -> (uint64 row * 1UL<Vector.index>, v))
+
+    Vector.fromCoordinateList (Vector.CoordinateList(uint64 rows * 1UL<Vector.dataLength>, vectorData))
+
+let reduceCols (op: 'a option -> 'a option -> 'a option) (matrix: SparseMatrix<'a>) : Vector.SparseVector<'a> =
+    let cols = matrix.ncols
+
+    let rec inner (size: uint64<storageSize>) (row: uint64<rowindex>) (col: uint64<colindex>) tree acc =
+        match tree with
+        | Node(nw, ne, sw, se) ->
+            let half = size / 2UL
+            let halfRow = (uint64 half) * 1UL<rowindex>
+            let halfCol = (uint64 half) * 1UL<colindex>
+            let acc = inner half row col nw acc
+            let acc = inner half row (col + halfCol) ne acc
+            let acc = inner half (row + halfRow) col sw acc
+            let acc = inner half (row + halfRow) (col + halfCol) se acc
+            acc
+        | Leaf(Dummy) -> acc
+        | Leaf(UserValue(None)) -> acc
+        | Leaf(UserValue(Some v)) -> (col, v) :: acc
+
+    let pairs =
+        inner matrix.storage.size 0UL<rowindex> 0UL<colindex> matrix.storage.data []
+
+    let grouped =
+        List.sortBy fst pairs
+        |> List.groupBy fst
+        |> List.map (fun (col, colSet) -> col, colSet |> List.map snd |> List.map Some)
+
+    let reduced =
+        grouped
+        |> List.map (fun (col, values) ->
+            match values with
+            | [] -> col, None
+            | head :: tail -> col, List.fold op head tail)
+
+    let data =
+        reduced
+        |> List.choose (fun (col, v) ->
+            match v with
+            | None -> None
+            | Some x -> Some(0UL<rowindex>, col, x))
+        |> List.sort
+
+    let vectorData =
+        data |> List.map (fun (_, col, v) -> (uint64 col * 1UL<Vector.index>, v))
+
+    Vector.fromCoordinateList (Vector.CoordinateList(uint64 cols * 1UL<Vector.dataLength>, vectorData))
+
+let kroneckerProduct
+    (matrix1: SparseMatrix<'a>)
+    (matrix2: SparseMatrix<'b>)
+    (f: 'a -> 'b -> 'c option)
+    : Result<SparseMatrix<'c>, string> =
+    let newRows = uint64 matrix1.nrows * uint64 matrix2.nrows * 1UL<nrows>
+    let newCols = uint64 matrix1.ncols * uint64 matrix2.ncols * 1UL<ncols>
+
+    let newSize =
+        getNearestUpperPowerOfTwo (max (uint64 newRows) (uint64 newCols))
+        * 1UL<storageSize>
+
+    let rec empty (size: uint64<storageSize>) =
+        match size with
+        | 1UL<storageSize> -> Leaf Dummy
+        | _ ->
+            let half = size / 2UL
+            let subtree = empty half
+            Node(subtree, subtree, subtree, subtree)
+
+    let rec insert (size: uint64<storageSize>) (row: uint64<rowindex>) (col: uint64<colindex>) value tree =
+        match size with
+        | 1UL<storageSize> -> Leaf(UserValue(value))
+        | _ ->
+            let half = size / 2UL
+            let halfRow = (uint64 half) * 1UL<rowindex>
+            let halfCol = (uint64 half) * 1UL<colindex>
+
+            match tree with
+            | Node(nw, ne, sw, se) ->
+                match (row < halfRow, col < halfCol) with
+                | true, true -> Node(insert half row col value nw, ne, sw, se)
+                | true, false -> Node(nw, insert half row (col - halfCol) value ne, sw, se)
+                | false, true -> Node(nw, ne, insert half (row - halfRow) col value sw, se)
+                | false, false -> Node(nw, ne, sw, insert half (row - halfRow) (col - halfCol) value se)
+            | _ ->
+                let emptySub = empty half
+
+                match (row < halfRow, col < halfCol) with
+                | true, true -> Node(insert half row col value emptySub, emptySub, emptySub, emptySub)
+                | true, false -> Node(emptySub, insert half row (col - halfCol) value emptySub, emptySub, emptySub)
+                | false, true -> Node(emptySub, emptySub, insert half (row - halfRow) col value emptySub, emptySub)
+                | false, false ->
+                    Node(emptySub, emptySub, emptySub, insert half (row - halfRow) (col - halfCol) value emptySub)
+
+    let foldQuadtree folder state size tree =
+        let rec inner rowOffset colOffset currentSize subTree acc =
+            match subTree with
+            | Leaf Dummy -> acc
+            | Leaf(UserValue None) -> acc
+            | Leaf(UserValue(Some value)) ->
+                let rec loop currRow currCol currentAcc =
+                    if currRow = rowOffset + currentSize then
+                        currentAcc
+                    elif currCol = colOffset + currentSize then
+                        loop (currRow + 1UL) colOffset currentAcc
+                    else
+                        loop currRow (currCol + 1UL) (folder currentAcc currRow currCol value)
+
+                loop rowOffset colOffset acc
+            | Node(nw, ne, sw, se) ->
+                let half = currentSize / 2UL
+                let acc1 = inner rowOffset colOffset half nw acc
+                let acc2 = inner rowOffset (colOffset + half) half ne acc1
+                let acc3 = inner (rowOffset + half) colOffset half sw acc2
+                inner (rowOffset + half) (colOffset + half) half se acc3
+
+        inner 0UL 0UL (uint64 size) tree state
+
+    let mat2Rows = uint64 matrix2.nrows
+    let mat2Cols = uint64 matrix2.ncols
+    let initialTree = empty newSize
+
+    let foldMatrixB acc rowMat1 colMat1 valMat1 =
+        foldQuadtree
+            (fun currentAcc rowMat2 colMat2 valMat2 ->
+                match f valMat1 valMat2 with
+                | Some computedVal ->
+                    let destRow = rowMat1 * mat2Rows + rowMat2
+                    let destCol = colMat1 * mat2Cols + colMat2
+                    insert newSize (destRow * 1UL<rowindex>) (destCol * 1UL<colindex>) (Some computedVal) currentAcc
+                | None -> currentAcc)
+            acc
+            matrix2.storage.size
+            matrix2.storage.data
+
+    let finalTree =
+        foldQuadtree
+            (fun acc rowMat1 colMat1 valMat1 -> foldMatrixB acc rowMat1 colMat1 valMat1)
+            initialTree
+            matrix1.storage.size
+            matrix1.storage.data
+
+    let rec count (size: uint64<storageSize>) tree =
+        match tree with
+        | Node(nw, ne, sw, se) ->
+            let half = size / 2UL
+            count half nw + count half ne + count half sw + count half se
+        | Leaf(UserValue(Some _)) -> 1UL<nvals>
+        | _ -> 0UL<nvals>
+
+    let nvals = count newSize finalTree
+
+    Ok(SparseMatrix(newRows, newCols, nvals, Storage(newSize, finalTree)))
diff --git a/QuadTree/Vector.fs b/QuadTree/Vector.fs
index b7bcb44..42961f2 100644
--- a/QuadTree/Vector.fs
+++ b/QuadTree/Vector.fs
@@ -307,6 +307,16 @@ let update (vector: SparseVector<_>) i v op =
         Error Error.InconsistentSizeOfArguments
 
 let fromCoordinateList (lst: CoordinateList<'a>) : SparseVector<'a> =
+    let unique =
+        lst.data
+        |> List.groupBy fst
+        |> List.map (fun (idx, entries) ->
+            let value = entries |> List.map snd |> List.last
+            (idx, value))
+
+    if unique |> List.exists (fun (idx, _) -> uint64 idx >= uint64 lst.length) then
+        failwith "Index out of range"
+
     let length = lst.length
     let nvals = (uint64 <| List.length lst.data) * 1UL<nvals>
     let storageSize = (getNearestUpperPowerOfTwo <| uint64 length) * 1UL<storageSize>
@@ -325,7 +335,7 @@ let fromCoordinateList (lst: CoordinateList<'a>) : SparseVector<'a> =
 
             mkNode left right, rCoordinates
 
-    let sortedCoordinates = List.sort lst.data
+    let sortedCoordinates = List.sort unique
 
     let tree, _ =
         traverse sortedCoordinates 0UL<index> ((uint64 storageSize) * 1UL<index>)
@@ -561,3 +571,91 @@ let scatter
                 | Error x -> Error x)
             (Ok w)
     | Error x -> Error Error.InconsistentStructureOfStorages
+
+
+let slice (_start: int) (_end: int) (vector: SparseVector<'a>) : Result<SparseVector<'a>, string> =
+    if _start < 0 then
+        Error "Start should be >= 0"
+    elif _end < 0 then
+        Error "End should be >= 0"
+    elif _start > int vector.length - 1 then
+        Error "Start is out of Vector length"
+    elif _end > int vector.length - 1 then
+        Error "End is out of Vector length"
+    elif _start > _end then
+        Error "End should be >= Start"
+    else
+        let startIdx = uint64 _start * 1UL<index>
+        let endIdx = uint64 _end * 1UL<index>
+        let newLength = uint64 (_end - _start + 1) * 1UL<dataLength>
+        let newSize = getNearestUpperPowerOfTwo (uint64 newLength) * 1UL<storageSize>
+
+        let rec cut (size: uint64<storageSize>) (pos: uint64<index>) tree =
+            let sizeIdx = (uint64 size) * 1UL<index>
+
+            match tree with
+            | Node(l, r) ->
+                let half = size / 2UL
+                let halfIdx = (uint64 half) * 1UL<index>
+                Node(cut half pos l, cut half (pos + halfIdx) r)
+            | Leaf(Dummy) -> Leaf Dummy
+            | Leaf(UserValue(v)) when pos >= startIdx && pos + sizeIdx - 1UL<index> <= endIdx -> Leaf(UserValue(v))
+            | _ -> Leaf Dummy
+
+        let cutTree = cut vector.storage.size 0UL<index> vector.storage.data
+
+        let rec empty (size: uint64<storageSize>) =
+            match size with
+            | 1UL<storageSize> -> Leaf Dummy
+            | _ ->
+                let half = size / 2UL
+                let subtree = empty half
+                Node(subtree, subtree)
+
+        let rec insert (size: uint64<storageSize>) (idx: uint64<index>) value tree =
+            match size with
+            | 1UL<storageSize> -> Leaf(UserValue(value))
+            | _ ->
+                let half = size / 2UL
+                let border = (uint64 half) * 1UL<index>
+
+                match tree with
+                | Node(left, right) ->
+                    if idx < border then
+                        Node(insert half idx value left, right)
+                    else
+                        Node(left, insert half (idx - border) value right)
+                | _ ->
+                    let emptySub = empty half
+
+                    if idx < border then
+                        Node(insert half idx value emptySub, emptySub)
+                    else
+                        Node(emptySub, insert half (idx - border) value emptySub)
+
+        let rec rebuild (size: uint64<storageSize>) (pos: uint64<index>) tree acc =
+            match tree with
+            | Leaf(Dummy) -> acc
+            | Leaf(UserValue(v)) ->
+                let newIdx = pos - startIdx
+                insert newSize newIdx v acc
+            | Node(left, right) ->
+                let half = size / 2UL
+                let border = (uint64 half) * 1UL<index>
+                let acc = rebuild half pos left acc
+                rebuild half (pos + border) right acc
+
+        let emptyTree = empty newSize
+        let shiftedTree = rebuild vector.storage.size 0UL<index> cutTree emptyTree
+
+        let rec count (size: uint64<storageSize>) tree =
+            match tree with
+            | Node(l, r) ->
+                let half = size / 2UL
+                count half l + count half r
+            | Leaf(UserValue(_)) -> (uint64 size) * 1UL<nvals>
+            | _ -> 0UL<nvals>
+
+        let nvals = count newSize shiftedTree
+
+        Ok(SparseVector(newLength, nvals, Storage(newSize, shiftedTree)))