diff --git a/docs/make.jl b/docs/make.jl
index 2417d00d..507d1db2 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -19,7 +19,8 @@ makedocs(
             "selectops.md"#,
             #"udfs.md"
         ],
-        "Utilities" => "utilities.md"
+        "Utilities" => "utilities.md",
+        "Examples" => "examples.md"
     ]
 )
 
diff --git a/docs/src/arrays.md b/docs/src/arrays.md
index 5ec56be8..614f54ac 100644
--- a/docs/src/arrays.md
+++ b/docs/src/arrays.md
@@ -1,5 +1,10 @@
 # Array Types
 
+```@meta
+DocTestSetup = :(using SuiteSparseGraphBLAS, SparseArrays)
+DocTestFilters = [r"memory: \d+ bytes"]
+```
+
 There are two primary array types in SuiteSparseGraphBLAS.jl: [`GBVector`](@ref) and [`GBMatrix`](@ref), as well as a few specialized versions of those array types. The full type hierarchy is:
 
 ```
@@ -33,17 +38,38 @@ Users should rarely need to directly interact with the underlying storage format
 ### Construction
 
 There are several methods to construct GBArrays. Shown here are empty construction, conversion from a dense matrix and a sparse matrix, and coordinate form with uniform or *iso* coefficients. 
-```@setup mat
-using SuiteSparseGraphBLAS
-using SparseArrays
-```
-```@repl mat
-x = GBMatrix{Bool}(20_000_000, 50_000)
-x = GBMatrix([[1,2] [3,4]])
-x = GBMatrix(sprand(100, 100, 0.5); fill = 0.0)
-x = GBMatrix(
-    rand(1:50_000, 5000), rand(1:500_000, 5000), 1, 500_000, 500_000
-)
+```jldoctest mat
+julia> x = GBMatrix{Bool}(20_000_000, 50_000)
+20000000x50000 GraphBLAS bool matrix, hypersparse by row, ints: 32/32/32
+  no entries, memory: 256 bytes
+
+julia> x = GBMatrix([[1,2] [3,4]])
+2x2 GraphBLAS int64_t matrix, full by col
+  4 entries, memory: 272 bytes
+
+    (1,1)   1
+    (2,1)   2
+    (1,2)   3
+    (2,2)   4
+
+julia> x = GBMatrix(sparse([1, 2, 3], [3, 1, 2], [1.0, 2.0, 3.0]); fill = 0.0)
+3x3 GraphBLAS double matrix, bitmap by col
+  3 entries, memory: 321 bytes
+
+    (2,1)    2
+    (3,2)    3
+    (1,3)    1
+
+julia> x = GBMatrix(
+           [1, 20_000, 400_000], [7, 30_000, 500_000], 1, 500_000, 500_000
+       )
+500000x500000 GraphBLAS int64_t matrix, hypersparse by row, ints: 32/32/32
+  3 entries, memory: 288 bytes
+  iso value:   1
+
+    (1,7)   1
+    (20000,30000)   1
+    (400000,500000)   1
 ```
 
 ```@docs
@@ -56,10 +82,26 @@ A `GBVector` is the one-dimensional equivalent of the `GBMatrix`, and internally
 
 ### Construction 
 
-```@repl mat
-v = GBVector{ComplexF32}(100)
-v = GBMatrix(rand(ComplexF64, 3); fill = nothing)
-v = GBVector(sprand(Bool, 100_000_000, 0.001))
+```jldoctest mat
+julia> v = GBVector{ComplexF32}(100)
+100x1 GraphBLAS float complex matrix, sparse by col, ints: 32/32/32
+  no entries, memory: 248 bytes
+
+julia> v = GBMatrix(ComplexF64[1, 2im, 3 + 4im]; fill = nothing)
+3x1 GraphBLAS double complex matrix, full by col
+  3 entries, memory: 288 bytes
+
+    (1,1)    1 + 0i
+    (2,1)    0 + 2i
+    (3,1)    3 + 4i
+
+julia> v = GBVector(sparsevec([1, 50_000_000, 100_000_000], true, 100_000_000))
+100000000x1 GraphBLAS bool matrix, sparse by col, ints: 64/64/64
+  3 entries, memory: 283 bytes
+
+    (1,1)   1
+    (50000000,1)   1
+    (100000000,1)   1
 ```
 
 ```@docs
@@ -81,15 +123,87 @@ The usual AbstractArray and SparseArray indexing capabilities are available. Inc
 
     However, many functions outside of GraphBLAS will throw an error if they receive `nothing` from an indexing operation. To accomodate these functions the user may set the fill value for an `AbstractGBArray` on construction and with [`setfill`](@ref) and [`setfill!`](@ref).
 
-```@repl mat
-A = GBMatrix([1,1,2,2,3,4,4,5,6,7,7,7], [2,4,5,7,6,1,3,6,3,3,4,5], [1:12...])
-SparseMatrixCSC(A)
-A[4]
-A[1,2]
-A[[1,3,5,7], :]
-A[1:2:7, :]
-A[:,:]
-A[:, 5]
+```jldoctest mat
+julia> A = GBMatrix([1,1,2,2,3,4,4,5,6,7,7,7], [2,4,5,7,6,1,3,6,3,3,4,5], [1:12...])
+7x7 GraphBLAS int64_t matrix, bitmap by row
+  12 entries, memory: 681 bytes
+
+    (1,2)   1
+    (1,4)   2
+    (2,5)   3
+    (2,7)   4
+    (3,6)   5
+    (4,1)   6
+    (4,3)   7
+    (5,6)   8
+    (6,3)   9
+    (7,3)   10
+    (7,4)   11
+    (7,5)   12
+
+julia> SparseMatrixCSC(A)
+7×7 SparseMatrixCSC{Int64, Int64} with 12 stored entries:
+ ⋅  1   ⋅   2   ⋅  ⋅  ⋅
+ ⋅  ⋅   ⋅   ⋅   3  ⋅  4
+ ⋅  ⋅   ⋅   ⋅   ⋅  5  ⋅
+ 6  ⋅   7   ⋅   ⋅  ⋅  ⋅
+ ⋅  ⋅   ⋅   ⋅   ⋅  8  ⋅
+ ⋅  ⋅   9   ⋅   ⋅  ⋅  ⋅
+ ⋅  ⋅  10  11  12  ⋅  ⋅
+
+julia> A[4]
+6
+
+julia> A[1,2]
+1
+
+julia> A[[1,3,5,7], :]
+4x7 GraphBLAS int64_t matrix, bitmap by row
+  7 entries, memory: 492 bytes
+
+    (1,2)   1
+    (1,4)   2
+    (2,6)   5
+    (3,6)   8
+    (4,3)   10
+    (4,4)   11
+    (4,5)   12
+
+julia> A[1:2:7, :]
+4x7 GraphBLAS int64_t matrix, bitmap by row
+  7 entries, memory: 492 bytes
+
+    (1,2)   1
+    (1,4)   2
+    (2,6)   5
+    (3,6)   8
+    (4,3)   10
+    (4,4)   11
+    (4,5)   12
+
+julia> A[:,:]
+7x7 GraphBLAS int64_t matrix, bitmap by row
+  12 entries, memory: 681 bytes
+
+    (1,2)   1
+    (1,4)   2
+    (2,5)   3
+    (2,7)   4
+    (3,6)   5
+    (4,1)   6
+    (4,3)   7
+    (5,6)   8
+    (6,3)   9
+    (7,3)   10
+    (7,4)   11
+    (7,5)   12
+
+julia> A[:, 5]
+7x1 GraphBLAS int64_t matrix, bitmap by col
+  2 entries, memory: 304 bytes
+
+    (2,1)   3
+    (7,1)   12
 ```
 
 The functionality illustrated above extends to `GBVector` as well.
diff --git a/docs/src/binaryops.md b/docs/src/binaryops.md
index 759e2033..e13c6b64 100644
--- a/docs/src/binaryops.md
+++ b/docs/src/binaryops.md
@@ -1,5 +1,10 @@
 # Binary Operators
 
+```@meta
+DocTestSetup = :(using SuiteSparseGraphBLAS)
+DocTestFilters = [r"memory: \d+ bytes", r"Ptr{.*}\(0x[0-9a-f]+\)"]
+```
+
 Binary operators are defined on three domains $D_1 \times D_2 \rightarrow D_3$.
 However, the vast majority of binary operators are defined on a single domain.
 
@@ -7,29 +12,83 @@ However, the vast majority of binary operators are defined on a single domain.
 
 In almost all cases you should pass Julia functions, which will be mapped to built-in operators, or used to create a new user-defined operator.
 
-```@repl
-using SuiteSparseGraphBLAS
+```jldoctest
+julia> x = GBMatrix([[1,2] [3,4]])
+2x2 GraphBLAS int64_t matrix, full by col
+  4 entries, memory: 272 bytes
+
+    (1,1)   1
+    (2,1)   2
+    (1,2)   3
+    (2,2)   4
+
+julia> x .+ x
+2x2 GraphBLAS int64_t matrix, full by row
+  4 entries, memory: 272 bytes
+
+    (1,1)   2
+    (1,2)   6
+    (2,1)   4
+    (2,2)   8
+
+julia> eadd(x, x, +)
+2x2 GraphBLAS int64_t matrix, full by row
+  4 entries, memory: 272 bytes
+
+    (1,1)   2
+    (1,2)   6
+    (2,1)   4
+    (2,2)   8
+
+julia> x .^ x
+2x2 GraphBLAS int64_t matrix, full by row
+  4 entries, memory: 272 bytes
+
+    (1,1)   1
+    (1,2)   27
+    (2,1)   4
+    (2,2)   256
+
+julia> emul(x, x, ^)
+2x2 GraphBLAS int64_t matrix, full by row
+  4 entries, memory: 272 bytes
+
+    (1,1)   1
+    (1,2)   27
+    (2,1)   4
+    (2,2)   256
 
-x = GBMatrix([[1,2] [3,4]])
+julia> x2 = Float64.(x)
+2x2 GraphBLAS double matrix, full by col
+  4 entries, memory: 272 bytes
 
-x .+ x
-eadd(x, x, +)
+    (1,1)    1
+    (2,1)    2
+    (1,2)    3
+    (2,2)    4
 
-x .^ x
-emul(x, x, ^)
+julia> eadd!(x2, x, x, +; accum=/)
+2x2 GraphBLAS double matrix, full by col
+  4 entries, memory: 272 bytes
 
-x2 = Float64.(x)
-eadd!(x2, x, x, +; accum=/)
+    (1,1)    0.5
+    (2,1)    0.5
+    (1,2)    0.5
+    (2,2)    0.5
 ```
 
 Internally functions are lowered like this:
 
-```@repl
-using SuiteSparseGraphBLAS
+```jldoctest
+julia> typedop = binaryop(+, Int64, Int64)
+SuiteSparseGraphBLAS.TypedBinaryOperator{typeof(+), Int64, Int64, Int64}(true, true, "GrB_PLUS_INT64", Ptr{SuiteSparseGraphBLAS.LibGraphBLAS.GB_BinaryOp_opaque}(0x0000736051461e20), +, nothing)
 
-typedop = binaryop(+, Int64, Int64)
+julia> eadd(GBVector([1,2]), GBVector([3,4]), typedop)
+2x1 GraphBLAS int64_t matrix, full by col
+  2 entries, memory: 256 bytes
 
-eadd(GBVector([1,2]), GBVector([3,4]), typedop)
+    (1,1)   4
+    (2,1)   6
 ```
 
 ## Built-Ins
diff --git a/docs/src/examples.md b/docs/src/examples.md
new file mode 100644
index 00000000..3b1798e1
--- /dev/null
+++ b/docs/src/examples.md
@@ -0,0 +1,71 @@
+# Examples
+
+```@meta
+DocTestSetup = quote
+    using SuiteSparseGraphBLAS
+    for f in ("Cohen", "Sandia", "BFSParents", "BellmanFord")
+        include(joinpath(pkgdir(SuiteSparseGraphBLAS), "examples", "$f.jl"))
+    end
+end
+DocTestFilters = [r"memory: \d+ bytes"]
+```
+
+The [`examples`](https://github.com/JuliaSparse/SuiteSparseGraphBLAS.jl/tree/v11/examples) directory of the repository contains graph algorithms written with SuiteSparseGraphBLAS.jl.
+They are not part of the package API: `include` a file to use the functions it defines.
+Each file links to a paper or reference describing the algorithm, and all of them are tested in `test/examples.jl`.
+
+Unless noted otherwise `A` is the adjacency matrix of an undirected graph with `n` vertices.
+
+| Algorithm | File | Function | Notes |
+|:----------|:-----|:---------|:------|
+| Triangle counting (Cohen) | `Cohen.jl` | `cohen(A)` | |
+| Triangle counting (Sandia) | `Sandia.jl` | `sandia(A)` | |
+| Triangle counting (CMU) | `CMU.jl` | `cmu(A, n)` | |
+| Triangles per vertex | `NodeWiseTriangleCount.jl` | `nodewise_tc(A)` | |
+| k-truss | `KTruss.jl` | `ktruss(A, k)` | Returns the truss with edge supports as values. |
+| BFS parents | `BFSParents.jl` | `bfs_parents(A, s, n)` | The parent of the source is `0`. |
+| Multi-source BFS parents | `BFSParentsMultiSource.jl` | `ms_bfsp(A, S, n, s)` | Row `i` of the `s × n` matrix `S` marks source `i`. |
+| Multi-source BFS levels | `BFSLevelsMultiSource.jl` | `ms_bfsl(A, S, n, s)` | `A` and `S` are `Bool`. Sources are level `1`, unreachable vertices `0`. |
+| Bellman-Ford | `BellmanFord.jl` | `bellmanford(A, s, n)` | `A[i,j]` is the weight of edge `i → j`, with explicit zeros on the diagonal. |
+| Floyd-Warshall | `FloydWarshall.jl` | `floydwarshall(A, n)` | Same input as Bellman-Ford. |
+| Prim's minimum spanning tree | `BasicPrim.jl`, `AdvancedPrim.jl` | `b_prim(A, n, zeros(n))`, `a_prim(A, n, zeros(n))` | `A` is dense, with `Inf` where there is no edge and zeros on the diagonal. |
+| Community detection by label propagation | `CDLP.jl` | `cdlp(A, n, t)` | Runs `t` iterations. Every vertex must have a neighbour. |
+
+The graph below has edges `1-2`, `1-3`, `2-3`, `2-4`, `3-4` and `4-5`, and two triangles, `{1,2,3}` and `{2,3,4}`.
+
+```jldoctest examples
+julia> src = [1, 1, 2, 2, 3, 4]; dst = [2, 3, 3, 4, 4, 5];
+
+julia> A = GBMatrix([src; dst], [dst; src], 1, 5, 5);
+
+julia> cohen(A), sandia(A)
+(2, 2)
+
+julia> bfs_parents(A, 1, 5)
+5x1 GraphBLAS int64_t matrix, bitmap by col
+  5 entries, memory: 288 bytes
+
+    (1,1)   0
+    (2,1)   1
+    (3,1)   1
+    (4,1)   2
+    (5,1)   4
+```
+
+With edge weights, and the explicit zero diagonal that Bellman-Ford expects:
+
+```jldoctest examples
+julia> w = [1.0, 4.0, 2.0, 7.0, 3.0, 1.0];
+
+julia> W = GBMatrix([src; dst; 1:5], [dst; src; 1:5], [w; w; zeros(5)]);
+
+julia> bellmanford(W, 1, 5)
+5x1 GraphBLAS double matrix, full by col
+  5 entries, memory: 280 bytes
+
+    (1,1)    0
+    (2,1)    1
+    (3,1)    3
+    (4,1)    6
+    (5,1)    7
+```
diff --git a/docs/src/index.md b/docs/src/index.md
index 59f04712..ca3425e4 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -1,5 +1,10 @@
 # SuiteSparseGraphBLAS.jl
 
+```@meta
+DocTestSetup = :(using SuiteSparseGraphBLAS)
+DocTestFilters = [r"memory: \d+ bytes"]
+```
+
 Fast sparse linear algebra is an essential part of the scientific computing toolkit. Outside of the usual applications, like differential equations, sparse linear algebra provides an elegant way to express graph algorithms on adjacency and incidence matrices. The GraphBLAS standard specifies a set of operations for computing sparse matrix graph algorithm in a vein similar to the BLAS or LAPACK standards.
 
 SuiteSparseGraphBLAS.jl is a blazing fast package for shared memory sparse matrix operations which wraps Tim Davis' SuiteSparse:GraphBLAS. If you use this package in your research please see [Citing](@ref).
@@ -41,21 +46,21 @@ The core SuiteSparseGraphBLAS.jl array types are `GBVector` and `GBMatrix` which
 !!! note "GBArray"
     These docs will often refer to the `GBArray` type, which is the union of `AbstractGBVector`, `AbstractGBMatrix` and their lazy Transpose objects.
 
-```@setup intro
-using SuiteSparseGraphBLAS
-using SparseArrays
-```
-
-```@repl intro
-# create a size 13 empty sparse vector with Float64 elements.
-v = GBVector{Float64}(13) 
-
-# create a 1000 x 1000 empty sparse matrix with ComplexF64 elements.
-A = GBMatrix{ComplexF64}(1000, 1000)
-
-# Non-stored values are equal to the fill value of A, 
-# which is by default zero(eltype(A))
-A[1,5] == getfill(A)
+```jldoctest intro
+julia> # create a size 13 empty sparse vector with Float64 elements.
+       v = GBVector{Float64}(13)
+13x1 GraphBLAS double matrix, sparse by col, ints: 32/32/32
+  no entries, memory: 248 bytes
+
+julia> # create a 1000 x 1000 empty sparse matrix with ComplexF64 elements.
+       A = GBMatrix{ComplexF64}(1000, 1000)
+1000x1000 GraphBLAS double complex matrix, hypersparse by row, ints: 32/32/32
+  no entries, memory: 256 bytes
+
+julia> # Non-stored values are equal to the fill value of A,
+       # which is by default zero(eltype(A))
+       A[1,5] == getfill(A)
+true
 ```
 
 Here we can already see several differences compared to `SparseArrays.SparseMatrixCSC`.
@@ -78,19 +83,43 @@ As noted in the example above, `A` has a fill-value, which defaults to `zero(elt
 For normal linear algebra this is a good choice for fill-value, however GraphBLAS is intended for graphs.
 To handle this `SuiteSparseGraphBLAS.jl` support `missing` for the fill value (this is an active area of development, and a new fill value that acts like the identity for most operations is better suited to this task).
 
-```@repl intro
-A[1, 1] == zero(eltype(A))
+```jldoctest intro
+julia> A[1, 1] == zero(eltype(A))
+true
+
+julia> B = setfill(A, missing) # no-copy alias
+1000x1000 GraphBLAS double complex matrix, hypersparse by row, ints: 32/32/32
+  no entries, memory: 256 bytes
 
-B = setfill(A, missing) # no-copy alias
-B[1, 1]
+julia> B[1, 1]
+missing
 ```
 
 An empty matrix and vector won't do us much good, so let's see how to construct the matrix and vector from the graphic above. Both `A` and `v` below are constructed from coordinate format or COO.
 
-```@repl intro
-A = GBMatrix([1,1,2,2,3,4,4,5,6,7,7,7], [2,4,5,7,6,1,3,6,3,3,4,5], [1:12...])
-
-v = GBVector([4], [10], 7)
+```jldoctest intro
+julia> A = GBMatrix([1,1,2,2,3,4,4,5,6,7,7,7], [2,4,5,7,6,1,3,6,3,3,4,5], [1:12...])
+7x7 GraphBLAS int64_t matrix, bitmap by row
+  12 entries, memory: 681 bytes
+
+    (1,2)   1
+    (1,4)   2
+    (2,5)   3
+    (2,7)   4
+    (3,6)   5
+    (4,1)   6
+    (4,3)   7
+    (5,6)   8
+    (6,3)   9
+    (7,3)   10
+    (7,4)   11
+    (7,5)   12
+
+julia> v = GBVector([4], [10], 7)
+7x1 GraphBLAS int64_t matrix, bitmap by col
+  1 entry, memory: 304 bytes
+
+    (4,1)   10
 ```
 
 ## GraphBLAS Operations
@@ -131,10 +160,40 @@ These built-in operators are *fast*, and should be used where possible. However,
 
 SuiteSparseGraphBLAS.jl will take care of operators behind the scenes, and in most cases users should pass in normal functions like `+` and `sin`. For example:
 
-```@repl intro
-emul(A, A, ^) # elementwise exponent
-
-map(sin, A)
+```jldoctest intro
+julia> emul(A, A, ^) # elementwise exponent
+7x7 GraphBLAS int64_t matrix, bitmap by row
+  12 entries, memory: 681 bytes
+
+    (1,2)   1
+    (1,4)   4
+    (2,5)   27
+    (2,7)   256
+    (3,6)   3125
+    (4,1)   46656
+    (4,3)   823543
+    (5,6)   16777216
+    (6,3)   387420489
+    (7,3)   10000000000
+    (7,4)   285311670611
+    (7,5)   8916100448256
+
+julia> map(sin, A)
+7x7 GraphBLAS double matrix, bitmap by row
+  12 entries, memory: 681 bytes
+
+    (1,2)    0.841471
+    (1,4)    0.909297
+    (2,5)    0.14112
+    (2,7)    -0.756802
+    (3,6)    -0.958924
+    (4,1)    -0.279415
+    (4,3)    0.656987
+    (5,6)    0.989358
+    (6,3)    0.412118
+    (7,3)    -0.544021
+    (7,4)    -0.99999
+    (7,5)    -0.536573
 ```
 
 Broadcasting functionality is also supported, `A .^ A` will lower to `emul(A, A, ^)`, and `sin.(A)` will lower to `apply(sin, A)`.
@@ -144,10 +203,27 @@ Matrix multiplication, which accepts a semiring, can be called with either `*(ma
 
 We can also use functions that are not already built into SuiteSparseGraphBLAS.jl:
 
-```@repl intro
-M = GBMatrix([[1,2] [3,4]])
-increment(x) = x + 1
-map(increment, M)
+```jldoctest intro
+julia> M = GBMatrix([[1,2] [3,4]])
+2x2 GraphBLAS int64_t matrix, full by col
+  4 entries, memory: 272 bytes
+
+    (1,1)   1
+    (2,1)   2
+    (1,2)   3
+    (2,2)   4
+
+julia> increment(x) = x + 1
+increment (generic function with 1 method)
+
+julia> map(increment, M)
+2x2 GraphBLAS int64_t matrix, full by col
+  4 entries, memory: 272 bytes
+
+    (1,1)   2
+    (2,1)   3
+    (1,2)   4
+    (2,2)   5
 ```
 
 Unfortunately this has a couple problems. The first is that it's slow.\
@@ -166,22 +242,52 @@ The methods are drawn from the LAGraph [repo](https://github.com/GraphBLAS/LAGra
 Input `A` must be a square, symmetric matrix with any element type.
 We'll test it using the matrix from the GBArray section above, which has two triangles in its undirected form.
 
-```@repl intro
-using SuiteSparseGraphBLAS: pair
-function cohen(A)
-  U = select(triu, A)
-  L = select(tril, A)
-  return reduce(+, *(L, U, (+, pair); mask=A)) ÷ 2
-end
-
-function sandia(A)
-  L = select(tril, A)
-  return reduce(+, *(L, L, (+, pair); mask=L))
-end
-
-M = eadd(A, A', +) #Make undirected/symmetric
-cohen(M)
-sandia(M)
+```jldoctest intro
+julia> using SuiteSparseGraphBLAS: pair
+
+julia> function cohen(A)
+         U = select(triu, A)
+         L = select(tril, A)
+         return reduce(+, *(L, U, (+, pair); mask=A)) ÷ 2
+       end
+cohen (generic function with 1 method)
+
+julia> function sandia(A)
+         L = select(tril, A)
+         return reduce(+, *(L, L, (+, pair); mask=L))
+       end
+sandia (generic function with 1 method)
+
+julia> M = eadd(A, A', +) #Make undirected/symmetric
+7x7 GraphBLAS int64_t matrix, bitmap by row
+  20 entries, memory: 681 bytes
+
+    (1,2)   1
+    (1,4)   8
+    (2,1)   1
+    (2,5)   3
+    (2,7)   4
+    (3,4)   7
+    (3,6)   14
+    (3,7)   10
+    (4,1)   8
+    (4,3)   7
+    (4,7)   11
+    (5,2)   3
+    (5,6)   8
+    (5,7)   12
+    (6,3)   14
+    (6,5)   8
+    (7,2)   4
+    (7,3)   10
+    (7,4)   11
+    (7,5)   12
+
+julia> cohen(M)
+2
+
+julia> sandia(M)
+2
 ```
 
 # Citing
diff --git a/docs/src/unaryops.md b/docs/src/unaryops.md
index 9facc5fe..6361ecaa 100644
--- a/docs/src/unaryops.md
+++ b/docs/src/unaryops.md
@@ -1,27 +1,53 @@
 # Unary Operators
 
+```@meta
+DocTestSetup = :(using SuiteSparseGraphBLAS)
+DocTestFilters = [r"memory: \d+ bytes", r"Ptr{.*}\(0x[0-9a-f]+\)"]
+```
+
 `UnaryOp`s are fairly straightforward $z = f(x)$ and their meaning should be clear from the name in most cases. 
 
 `UnaryOp`s are used only in the `map` and [`apply`](@ref) functions.
 
-```@repl
-using SuiteSparseGraphBLAS
+```jldoctest
+julia> x = GBVector([1.5, 0, pi])
+3x1 GraphBLAS double matrix, full by col
+  3 entries, memory: 264 bytes
+
+    (1,1)    1.5
+    (2,1)    0
+    (3,1)    3.14159
+
+julia> y = map(sin, x)
+3x1 GraphBLAS double matrix, full by col
+  3 entries, memory: 264 bytes
 
-x = GBVector([1.5, 0, pi])
+    (1,1)    0.997495
+    (2,1)    0
+    (3,1)    1.22465e-16
 
-y = map(sin, x)
+julia> map(asin, y)
+3x1 GraphBLAS double matrix, full by col
+  3 entries, memory: 264 bytes
 
-map(asin, y)
+    (1,1)    1.5
+    (2,1)    0
+    (3,1)    1.22465e-16
 ```
 
 Internally functions are lowered like this:
 
-```@repl
-using SuiteSparseGraphBLAS
+```jldoctest
+julia> op = unaryop(sin, Float64)
+SuiteSparseGraphBLAS.TypedUnaryOperator{typeof(sin), Float64, Float64}(true, true, "GxB_SIN_FP64", Ptr{SuiteSparseGraphBLAS.LibGraphBLAS.GB_UnaryOp_opaque}(0x0000736051454500), sin, nothing)
 
-op = unaryop(sin, Float64)
+julia> map(op, GBVector([1.5, 0, pi]))
+3x1 GraphBLAS double matrix, full by col
+  3 entries, memory: 264 bytes
 
-map(op, GBVector([1.5, 0, pi]))
+    (1,1)    0.997495
+    (2,1)    0
+    (3,1)    1.22465e-16
 ```
 
 ## Built-Ins
diff --git a/examples/AdvancedPrim.jl b/examples/AdvancedPrim.jl
new file mode 100644
index 00000000..8001f5c1
--- /dev/null
+++ b/examples/AdvancedPrim.jl
@@ -0,0 +1,47 @@
+#=
+AdvancedPrim:
+- Julia version: 
+- Author: samuel
+- Date: 2021-09-10
+- Reference: https://en.wikipedia.org/wiki/Prim%27s_algorithm
+=#
+
+using SuiteSparseGraphBLAS
+using SparseArrays
+using LinearAlgebra
+
+function stoppingCondition(m, n)
+
+infinito = true
+
+for i = 1:n
+        if(m[i] == 0.0)
+            return false
+    end
+end
+    return infinito
+end
+
+#A is the input matrix (dense, Inf where there is no edge and zeros on the diagonal), n is the number of nodes,
+#m is a vector of n zeros (Inf marks nodes already in the tree). Output: minimum spanning tree cost and edge weights in the order they were added
+function a_prim(A, n, m)
+
+    d = GBVector{Float64}(n)
+
+    weight = 0.0
+
+    d = A[1,:]
+    mst = Float64[] #minimum spanning tree
+
+    while(stoppingCondition(m, n) == false)
+
+        u = argmin(m+d)
+        m[u] = Inf
+        push!(mst, d[u])
+        weight = weight + d[u]
+        d = emul(d, A[u,:], min)
+    end
+
+    return weight, mst
+
+end
\ No newline at end of file
diff --git a/examples/BFSLevelsMultiSource.jl b/examples/BFSLevelsMultiSource.jl
new file mode 100644
index 00000000..b44d07b4
--- /dev/null
+++ b/examples/BFSLevelsMultiSource.jl
@@ -0,0 +1,49 @@
+#=
+BFSLevelsMultiSource:
+- Julia version: 1.6.2
+- Author: samuel
+- Date: 2021-09-10
+- Reference: https://en.wikipedia.org/wiki/Breadth-first_search
+=#
+
+using SuiteSparseGraphBLAS
+using SparseArrays
+using LinearAlgebra
+using SuiteSparseGraphBLAS: Complement, ∨, ∧
+
+#A is the input matrix, S is the source matrix, and n is the number of nodes, s is the number of sources in the matrix S, Output: Levels matrix
+function ms_bfsl(A, S, n, s)
+
+    #this vector allows to check if a value is already present inside the distance matrix,
+    #if it's present we don't want it to be changed
+    levels = GBVector{Int64}(n)
+        for i = 1:n
+            levels[i] = i
+        end
+
+    #distance matrix
+    distance = GBMatrix{Int64}(s, n)
+            for i = 1:s
+                for j = 1:n
+                    distance[i, j] = 0
+                end
+            end
+
+    for level = 1:n-1
+
+            for i = 1:s
+                for j = 1:n
+                    if(S[i, j] !== nothing && S[i, j] != 0 && distance[i, j] ∉ levels) #*
+                    distance[i, j] = level #*
+                    end
+                end
+            end
+
+
+            S = *(S, A, (∨, ∧), mask=Complement(distance), desc=Descriptor(replace_output=true))
+
+    end
+
+    return distance
+
+end
\ No newline at end of file
diff --git a/examples/BFSParents.jl b/examples/BFSParents.jl
new file mode 100644
index 00000000..3a27d434
--- /dev/null
+++ b/examples/BFSParents.jl
@@ -0,0 +1,75 @@
+#=
+BFSParents:
+- Julia version: 1.6.2
+- Author: samuel
+- Date: 2021-09-10
+- Reference: https://en.wikipedia.org/wiki/Breadth-first_search
+=#
+
+using SuiteSparseGraphBLAS
+using SparseArrays
+using LinearAlgebra
+using SuiteSparseGraphBLAS: Structural, Complement, second
+
+#function to insert inside the result array
+function insert(p, result)
+    for i=1:length(p)
+        if(p[i]!=nothing && result[i]==nothing)
+            result[i] = p[i]
+        end
+    end
+
+    return result
+
+end
+
+#function to update visited nodes
+function updateVisited(f, visited)
+    for i=1:length(f)
+        if(f[i]!=nothing)
+            visited[i] = true
+        end
+    end
+
+    return visited
+
+end
+
+#A is the input matrix, s is the source node, and n is the number of nodes, Output: parent matrix
+function bfs_parents(A, s, n)
+
+    index = GBVector{Int64}(n; fill=nothing)
+        for i = 1:n
+            index[i] = i
+        end
+
+    #wavefront vector -> w[s]=source node insrted in wavefront
+    f = GBVector{Int64}(n; fill=nothing)
+    f[s] = s
+
+    #parent vector -> p[s]=1 because the source is already visited
+    p = GBVector{Int64}(n; fill=nothing)
+
+    #result vector -> result[s]=0 because the parent of the source is considered node 0
+    result = GBVector{Int64}(n; fill=nothing)
+    result[s] = 0
+
+    #visited vector -> v[s]=1 because the source is already visited
+    visited = GBVector{Bool}(n; fill=nothing)
+    visited[s] = true
+
+    for i = 1:n-1
+            f = *(A', f, (min, second), mask=Complement(Structural(p)))
+            p = copy(f)
+            f = mask(index, Structural(f))
+            result = insert(p, result)
+            visited = updateVisited(f, visited)
+            if(nnz(visited) == n)
+                break
+            end
+
+    end
+
+    return result
+
+end
\ No newline at end of file
diff --git a/examples/BFSParentsMultiSource.jl b/examples/BFSParentsMultiSource.jl
new file mode 100644
index 00000000..706e1d21
--- /dev/null
+++ b/examples/BFSParentsMultiSource.jl
@@ -0,0 +1,68 @@
+#=
+BFSParentsMultiSource:
+- Julia version: 1.6.2
+- Author: samuel
+- Date: 2021-09-10
+- Reference: https://en.wikipedia.org/wiki/Breadth-first_search
+=#
+
+using SuiteSparseGraphBLAS
+using SparseArrays
+using LinearAlgebra
+using SuiteSparseGraphBLAS: Structural, Complement
+
+function insert(P, R, s, n)
+    for i = 1:s
+        for j = 1:n
+            if(P[i, j] !== nothing && R[i, j] === nothing) #*
+                R[i, j] = P[i, j]
+            end
+        end
+    end
+
+    return R
+
+end
+
+#A is the input matrix, S are the sources, n is the number of nodes, s is the number of sources, Output: Parents matrix
+function ms_bfsp(A, S, n, s)
+
+    index = GBMatrix{Int64}(s, n; fill=nothing)
+            for i = 1:n
+                for j = 1:s
+                    index[j, i] = i
+                end
+            end
+
+    #frontier matrix -> F=S source nodes insrted in frontier
+    F = mask(index, Structural(S))
+
+    #parent matrix -> P[S]=S because the source is already visited
+    P = F
+
+    #result matrix -> R[S]=0 because the parents of the source are considered node 0
+    R = GBMatrix{Int64}(s, n; fill=nothing)
+    LR = GBMatrix{Int64}(s, n; fill=nothing)
+    for i = 1:n
+        for j = 1:s
+            if(S[j, i] !== nothing && S[j, i] != 0)
+                R[j, i] = 0
+            end
+        end
+    end
+    for _ ∈ 1:n-1
+            F = *(F, A, (min, first), mask=Complement(P))
+            P = copy(F)
+            F = mask(index, Structural(F))
+            R = insert(P, R, s, n)
+            if(R == LR)
+                break
+            end
+            LR = copy(R)
+
+    end
+
+    return R
+
+end
+
diff --git a/examples/BasicPrim.jl b/examples/BasicPrim.jl
new file mode 100644
index 00000000..65ce2f01
--- /dev/null
+++ b/examples/BasicPrim.jl
@@ -0,0 +1,43 @@
+#=
+BasicPrim:
+- Julia version: 1.6.2
+- Author: samuel
+- Date: 2021-09-10
+- Reference: https://en.wikipedia.org/wiki/Prim%27s_algorithm
+=#
+
+using SuiteSparseGraphBLAS
+using SparseArrays
+
+function stoppingCondition(m, n)
+
+infinito = true
+
+for i = 1:n
+        if(m[i] == 0.0)
+            return false
+    end
+end
+    return infinito
+end
+
+#A is the input matrix (dense, Inf where there is no edge and zeros on the diagonal), n is the number of nodes,
+#m is a vector of n zeros (Inf marks nodes already in the tree) Output: minimum spanning tree cost
+function b_prim(A, n, m)
+
+    d = GBVector{Float64}(n)
+
+    weight = 0.0
+
+    d = A[1,:]
+    while(stoppingCondition(m, n) == false)
+
+        u = argmin(m+d)
+        m[u] = Inf
+        weight = weight + d[u]
+        d = min.(d, A[u,:])
+    end
+
+    return weight
+
+end
\ No newline at end of file
diff --git a/examples/BellmanFord.jl b/examples/BellmanFord.jl
new file mode 100644
index 00000000..8f56c50f
--- /dev/null
+++ b/examples/BellmanFord.jl
@@ -0,0 +1,29 @@
+#=
+BellmanFord:
+- Julia version: 1.6.2
+- Author: samuel
+- Date: 2021-09-10
+- Reference: https://en.wikipedia.org/wiki/Bellman%E2%80%93Ford_algorithm
+=#
+
+using SuiteSparseGraphBLAS
+using SparseArrays
+
+#A is the input matrix (A[i,j] is the weight of edge i->j, with explicit zeros on the diagonal),
+#s is the source node, and n is the number of nodes, Output: Shortest path from source node
+function bellmanford(A, s, n)
+
+    #vector init +inf
+    d = GBVector{Float64}(n)
+        for i = 1:n
+            d[i] = Inf
+        end
+
+    d[s]=0.0
+    for _ ∈ 1:n
+           d = *(A', d, (min, +))
+    end
+
+    return d
+
+end
\ No newline at end of file
diff --git a/examples/CDLP.jl b/examples/CDLP.jl
new file mode 100644
index 00000000..788a08fa
--- /dev/null
+++ b/examples/CDLP.jl
@@ -0,0 +1,64 @@
+#=
+CDLP:
+- Julia version: 1.6.2
+- Author: samuel
+- Date: 2021-09-10
+- Reference: https://doi.org/10.1103/PhysRevE.76.036106 (CDLP as defined in LDBC Graphalytics, https://arxiv.org/abs/2011.15028)
+=#
+
+using SuiteSparseGraphBLAS
+using SparseArrays
+using LinearAlgebra
+
+#minMode function: given a collection [1, 5, 5, 3, 3, 7] will return 3 because the highest repetition is 2
+#and the lowest most repeated number is 3
+function minMode(c)
+
+    temp = filter(x -> x!=0, c)
+    count = Dict{eltype(temp), Int}()
+    for x in temp
+        count[x] = get(count, x, 0) + 1
+    end
+    repetitions = maximum(values(count))
+    return minimum(key for (key, value) in count if value == repetitions)
+end
+
+#can be useful to tweak the algorithm
+function diag_conversion(lab, n)
+for i = 1:n
+        for j = 1:n
+            if(i==j)
+                lab[i,j]=i
+            else
+                lab[i,j]=0
+            end
+        end
+    end
+    return lab
+end
+
+#A is the input matrix, n is the number of nodes, t the number of iterations Output: label propagation matrix
+function cdlp(A, n, t)
+
+    lab = GBMatrix{Int64}(n,n)
+
+    for i = 1:n
+        for j = 1:n
+            if(i==j)
+                lab[i,j]=i
+            else
+                lab[i,j]=0
+            end
+        end
+    end
+
+    for k = 1:t
+            F = A * lab
+            for i = 1:n
+            r = nonzeros(F[i,:])
+            lab[i, i] = minMode(r)
+        end
+    end
+    return lab
+
+end
\ No newline at end of file
diff --git a/examples/CMU.jl b/examples/CMU.jl
new file mode 100644
index 00000000..cb8dec54
--- /dev/null
+++ b/examples/CMU.jl
@@ -0,0 +1,27 @@
+#=
+CMU:
+- Julia version: 1.6.2
+- Author: samuel
+- Date: 2021-09-10
+- Reference: https://doi.org/10.1109/HPEC.2017.8091046
+=#
+
+using SuiteSparseGraphBLAS
+using SparseArrays
+using LinearAlgebra
+
+#A is the input matrix, n is the number of nodes in the original matrix, Output: Triangle Count
+function cmu(A, n)
+    t = 0
+    for i ∈ 2:n-1
+        #rpva
+        A20 = A[i+1:n,begin:i]
+        a10 = A[begin:i, i]
+        a12 = A[i, i+1:n]
+
+        partial1 = a10' * A20'
+        partial2 = partial1 * a12'
+        t += partial2[1]
+    end
+    return t
+end
\ No newline at end of file
diff --git a/examples/Cohen.jl b/examples/Cohen.jl
new file mode 100644
index 00000000..c692dd0d
--- /dev/null
+++ b/examples/Cohen.jl
@@ -0,0 +1,19 @@
+#=
+Cohen:
+- Julia version: 1.6.2
+- Author: samuel
+- Date: 2021-09-10
+- Reference: https://doi.org/10.1109/MCSE.2009.120, https://doi.org/10.1109/IPDPSW.2015.75
+=#
+
+using SuiteSparseGraphBLAS
+using SparseArrays
+using LinearAlgebra
+using SuiteSparseGraphBLAS: pair
+
+##A is the input matrix, Output: Triangle Count
+function cohen(A)
+    U = triu(A)
+    L = tril(A)
+    return reduce(+, *(L, U, (+, pair); mask=A)) ÷ 2
+end
\ No newline at end of file
diff --git a/examples/FloydWarshall.jl b/examples/FloydWarshall.jl
new file mode 100644
index 00000000..d6f6405c
--- /dev/null
+++ b/examples/FloydWarshall.jl
@@ -0,0 +1,25 @@
+#=
+FloydWarshall:
+- Julia version: 1.6.2
+- Author: samuel
+- Date: 2021-09-10
+- Reference: https://en.wikipedia.org/wiki/Floyd%E2%80%93Warshall_algorithm
+=#
+
+using SuiteSparseGraphBLAS
+using SparseArrays
+
+#A is the input matrix for Floyd-Warshall multiple sources multiple destinations
+#(A[i,j] is the weight of edge i->j, with explicit zeros on the diagonal), Output: shortest paths
+function floydwarshall(A, n)
+
+    #distance matrix = input matrix
+    D = A
+
+    for i = 1:n
+            partial = *(D[:, i:i], D[i:i, :], (min, +))
+            D = eadd(D, partial, min)
+    end
+    return D
+
+end
\ No newline at end of file
diff --git a/examples/KTruss.jl b/examples/KTruss.jl
new file mode 100644
index 00000000..f33b28a9
--- /dev/null
+++ b/examples/KTruss.jl
@@ -0,0 +1,25 @@
+#=
+KTruss:
+- Julia version: 1.6.2
+- Author: samuel
+- Date: 2021-09-10
+- Reference: https://doi.org/10.1109/HPEC.2018.8547538
+=#
+
+using SuiteSparseGraphBLAS
+using SparseArrays
+using LinearAlgebra
+using SuiteSparseGraphBLAS: ∧
+
+#A is the original matrix, k is the requested truss, Output: K-Truss graph (values are edge supports)
+function ktruss(A, k)
+
+    C = A
+    while true
+       nz = nnz(C)
+       C = *(C, C, (+, ∧), mask=C)
+       C = select(>=, C, k-2)
+       nnz(C) == nz && return C
+    end
+
+end
\ No newline at end of file
diff --git a/examples/NodeWiseTriangleCount.jl b/examples/NodeWiseTriangleCount.jl
new file mode 100644
index 00000000..68d905f3
--- /dev/null
+++ b/examples/NodeWiseTriangleCount.jl
@@ -0,0 +1,20 @@
+#=
+NodeWiseTriangleCount:
+- Julia version: 1.6.2
+- Author: samuel
+- Date: 2021-09-10
+- Reference: https://doi.org/10.1109/HPEC.2018.8547538
+=#
+
+using SuiteSparseGraphBLAS
+using SparseArrays
+using LinearAlgebra
+
+#A is the input matrix. Output: triangles originated by each node
+function nodewise_tc(A)
+
+    partial = *(A, A, (+, *), mask=A)
+    t = reduce(+, partial, dims=2)
+    return t .÷ 2
+
+end
\ No newline at end of file
diff --git a/examples/Sandia.jl b/examples/Sandia.jl
new file mode 100644
index 00000000..0daeabdf
--- /dev/null
+++ b/examples/Sandia.jl
@@ -0,0 +1,20 @@
+#=
+Sandia:
+- Julia version: 1.6.2
+- Author: samuel
+- Date: 2021-09-10
+- Reference: https://doi.org/10.1109/HPEC.2017.8091043
+=#
+
+using SuiteSparseGraphBLAS
+using SparseArrays
+using LinearAlgebra
+
+#A is the input matrix, Output: Triangle Count
+function sandia(A)
+
+    L = tril(A)
+    C = *(L, L, (+, *); mask=L)
+    return reduce(+, C, dims=:)
+
+end
\ No newline at end of file
diff --git a/test/examples.jl b/test/examples.jl
new file mode 100644
index 00000000..730aacc4
--- /dev/null
+++ b/test/examples.jl
@@ -0,0 +1,62 @@
+module ExamplesTest
+using SuiteSparseGraphBLAS, SparseArrays, LinearAlgebra, Test
+
+for f in ["AdvancedPrim", "BasicPrim", "BellmanFord", "BFSLevelsMultiSource", "BFSParents",
+          "BFSParentsMultiSource", "CDLP", "CMU", "Cohen", "FloydWarshall", "KTruss",
+          "NodeWiseTriangleCount", "Sandia"]
+    include(joinpath(@__DIR__, "..", "examples", "$f.jl"))
+end
+
+# Undirected graph with edges 1-2, 1-3, 2-3, 2-4, 3-4, 4-5 and triangles {1,2,3}, {2,3,4}.
+I = [1, 1, 2, 2, 3, 4]; J = [2, 3, 3, 4, 4, 5]
+sym(v) = sparse([I; J], [J; I], [v; v], 5, 5)
+A = GBMatrix(sym(ones(Int, 6)))
+w = [1.0, 4.0, 2.0, 7.0, 3.0, 1.0]
+# Dense weights: Inf where there is no edge, zeros on the diagonal.
+W = fill(Inf, 5, 5); W[diagind(W)] .= 0
+for (i, j, x) in zip(I, J, w); W[i, j] = W[j, i] = x; end
+
+@testset "triangle counting" begin
+    @test cohen(A) == sandia(A) == cmu(A, 5) == 2
+    @test Vector(nodewise_tc(A)) == [1, 2, 2, 1, 0]
+    # K4 on 1:4 is the 4-truss of K4 plus the pendant edge 4-5.
+    K = GBMatrix(sparse([1, 1, 1, 2, 2, 3, 4], [2, 3, 4, 3, 4, 4, 5], 1, 5, 5))
+    K = K + K'
+    @test nnz(ktruss(A, 3)) == 10 && ktruss(A, 3)[4, 5] == 0
+    @test nnz(ktruss(K, 4)) == 12 && nnz(ktruss(K, 5)) == 0
+end
+
+@testset "shortest paths" begin
+    D = copy(W)
+    for k in 1:5, i in 1:5, j in 1:5
+        D[i, j] = min(D[i, j], D[i, k] + D[k, j])
+    end
+    G = GBMatrix(sparse([I; J; 1:5], [J; I; 1:5], [w; w; zeros(5)]))
+    @test Matrix(floydwarshall(G, 5)) == D
+    @test Vector(bellmanford(G, 1, 5)) == D[1, :] == [0, 1, 3, 6, 7]
+    @test Vector(bellmanford(G, 5, 5)) == D[5, :]
+end
+
+@testset "Prim" begin
+    @test b_prim(GBMatrix(W), 5, zeros(5)) == 7
+    @test a_prim(GBMatrix(W), 5, zeros(5)) == (7, [0, 1, 2, 3, 1])
+end
+
+@testset "BFS" begin
+    @test collect(bfs_parents(A, 1, 5)) == [0, 1, 1, 2, 4]
+    @test collect(bfs_parents(A, 5, 5)) == [2, 4, 4, 5, 0]
+    S = GBMatrix(sparse([1, 2], [1, 5], 1, 2, 5))
+    @test collect(ms_bfsp(A, S, 5, 2)) == [0 1 1 2 4; 2 4 4 5 0]
+    B = GBMatrix(sym(trues(6)))
+    @test Matrix(ms_bfsl(B, GBMatrix(sparse([1, 2], [1, 5], true, 2, 5)), 5, 2)) ==
+        [1 2 2 3 4; 4 3 3 2 1]
+end
+
+@testset "CDLP" begin
+    # Two disjoint triangles converge to the smallest label in each.
+    T = GBMatrix(sparse([1, 1, 2, 4, 4, 5], [2, 3, 3, 5, 6, 6], 1, 6, 6))
+    T = T + T'
+    @test diag(Matrix(cdlp(T, 6, 2))) == [1, 1, 1, 4, 4, 4]
+    @test diag(Matrix(cdlp(A, 5, 1))) == [2, 1, 1, 2, 4]
+end
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 3b044e20..4d50a992 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -85,4 +85,5 @@ println("$(SuiteSparseGraphBLAS.get_lib())")
     include_test("solvers/umfpack.jl")
     include_test("solvers/cholmod.jl")
     include_test("serialization.jl")
+    include_test("examples.jl")
 end