sufficient coverage of cholmod, wrappers still very broken. improve broadcasting slightly

Author: rayegun

src/SuiteSparseGraphBLAS.jl

@@ -101,6 +101,7 @@ include("chainrules/constructorrules.jl")
 include("serialization.jl")
 
 #EXPERIMENTAL
+include("linalg.jl")
 include("misc.jl")
 include("mmread.jl")
 include("iterator.jl")
@@ -133,7 +134,7 @@ export extract, extract!, subassign!, assign!, hvcat! #array functions
 
 #operations
 export select, select!, eadd, eadd!, emul, emul!, gbtranspose, gbtranspose!,
-gbrand, eunion, eunion!, mask, mask!, apply, apply!, setfill, setfill!
+gbrand, eunion, eunion!, mask, mask!, apply, apply!, setfill, setfill!, gbrandn
 # Reexports from LinAlg
 export diag, diagm, mul!, kron, kron!, transpose, reduce, tril, triu
 

src/abstractgbarray.jl

@@ -247,12 +247,12 @@ function build!(A::AbstractGBMatrix{T}, I::AbstractVector, J::AbstractVector, x:
     @wraperror LibGraphBLAS.GxB_Matrix_build_Scalar(
         A,
         Vector{LibGraphBLAS.GrB_Index}(decrement!(I)),
-        Vector{LibGraphBLAS.GrB_Index}(decrement!(J)),
+        Vector{LibGraphBLAS.GrB_Index}(I !== J ? decrement!(J) : J),
         x,
         length(I)
     )
     increment!(I)
-    increment!(J)
+    I !== J && (increment!(J))
     return A
 end
 
@@ -296,7 +296,7 @@ for T ∈ valid_vec
             length(X) == length(I) == length(J) ||
                 DimensionMismatch("I, J and X must have the same length")
             decrement!(I)
-            decrement!(J)
+            I !== J && (decrement!(J))
             @wraperror LibGraphBLAS.$func(
                 A,
                 I,
@@ -306,7 +306,7 @@ for T ∈ valid_vec
                 combine
             )
             increment!(I)
-            increment!(J)
+            I !== J && (increment!(J))
             return A
         end
     end
@@ -535,7 +535,7 @@ function subassign!(
     desc = _handledescriptor(desc; out=C, in1=A)
     mask = _handlemask!(desc, mask)
     I = decrement!(I)
-    J = decrement!(J)
+    I !== J && (J = decrement!(J))
     rereshape = false
     sz1 = size(A, 1)
     if !(eltype(A) <: valid_union) || !(eltype(C) <: valid_union)
@@ -553,7 +553,7 @@ function subassign!(
         parent(A), true, sz1, 1, C_NULL)
     end
     increment!(I)
-    increment!(J)
+    I !== J && (increment!(J))
     return A
 end
 
@@ -565,12 +565,12 @@ function subassign!(C::AbstractGBArray{T}, x, I, J;
     I, ni = idx(I)
     J, nj = idx(J)
     I = decrement!(I)
-    J = decrement!(J)
+    I !== J && (J = decrement!(J))
     desc = _handledescriptor(desc; out=C)
     mask = _handlemask!(desc, mask)
     _subassign(C, x, I, ni, J, nj, mask, _handleaccum(accum, storedeltype(C)), desc)
     increment!(I)
-    increment!(J)
+    I !== J && (decrement!(J))
     return x
 end
 
@@ -638,13 +638,13 @@ function assign!(
     desc = _handledescriptor(desc; in1=A, out=C)
     mask = _handlemask!(desc, mask)
     I = decrement!(I)
-    J = decrement!(J)
+    I !== J && (J = decrement!(J))
     if !(eltype(A) <: valid_union) || !(eltype(C) <: valid_union)
         A = LinearAlgebra.copy_oftype(A, eltype(C))
     end
     @wraperror LibGraphBLAS.GrB_Matrix_assign(C, mask, _handleaccum(accum, storedeltype(C)), parent(A), I, ni, J, nj, desc)
     increment!(I)
-    increment!(J)
+    I !== J && (decrement!(J))
     return A
 end
 
@@ -656,12 +656,12 @@ function assign!(C::AbstractGBArray{T}, x, I, J;
     I, ni = idx(I)
     J, nj = idx(J)
     I = decrement!(I)
-    J = decrement!(J)
+    I !== J && (J = decrement!(J))
     desc = _handledescriptor(desc; out=C)
     mask = _handlemask!(desc, mask)
     _assign(C, x, I, ni, J, nj, mask, _handleaccum(accum, storedeltype(C)), desc)
     increment!(I)
-    increment!(J)
+    I !== J && (decrement!(J))
     return x
 end
 
@@ -719,7 +719,7 @@ function GBDiagonal!(C::AbstractGBMatrix, v::AbstractGBVector, k::Integer=0; des
     return C
 end
 function GBDiagonal!(C::AbstractGBMatrix{T}, v::AbstractVector, k::Integer=0; desc = nothing) where T
-    v2 = GBShallowVector(convert(DenseVector{T}, v))
+    v2 = GBShallowVector(convert(Vector{T}, v))
     GBDiagonal!(C, v2, k; desc)
 end
 function GBDiagonal!(C::AbstractGBMatrix, D::Diagonal; desc = nothing)
@@ -736,7 +736,6 @@ function GBDiagonal(v::AbstractGBVector, k::Integer=0; desc = nothing)
     GBDiagonal!(C, v, k; desc)
 end
 
-
 # Type dependent functions build, setindex, getindex, and findnz:
 for T ∈ valid_vec
     if T ∈ gxb_vec

src/iterator.jl

@@ -1,9 +1,16 @@
 # TODO: Make this file less loathsome.
 
+# These iterators are only somewhat useful.
+# The biggest issue is that we can't, at the moment, use one to change values in another.
+# So in place map! is doable like this, as is just simply getting values.
+# But if you wanted to say `C[I] = <transform>(A[I])` even with identical structures
+# there's no method. Non-identical patterns are off the table entirely.
+
 abstract type IndexIteratorType end
 struct IndicesIterator <: IndexIteratorType end # return the indices as integers
 struct NeighborIterator <: IndexIteratorType end # Used only for Vector iterators
 struct NoIndexIterator <: IndexIteratorType end # Used when we don't want the indices, just the value
+struct IteratorIterator <: IndexIteratorType end # return the iterator object itself. Spooky!
 # just returns the free index, only useful for single column/row iterators.
 
 abstract type AbstractGBIterator end
@@ -143,13 +150,20 @@ end
     return unsafe_load(Ptr{T}(I.p.Ax[]), I.p.iso[] ? 1 : I.p.p[] + 1)
 end
 
+@inline function setval(I::GxBIterator{<:Any, T}, x::T) where T
+    I.p.iso[] && throw(ArgumentError("Cannot set value of iso valued matrix using iterator."))
+    unsafe_store!(Ptr{T}(I.p.Ax[]), x, I.p.p[] + 1)
+end
+
 # TODO: Inelegant
 @inline get_element(I::GxBIterator{<:Any, <:Any, true, IndicesIterator()}) = ((getrow(I), getcol(I)), getval(I))
 @inline get_element(I::GxBIterator{<:Any, <:Any, false, IndicesIterator()}) = (getrow(I), getcol(I))
 @inline get_element(I::GxBIterator{<:Any, <:Any, true, NeighborIterator()}) = (increment(_rc_geti(I)), getval(I))
 @inline get_element(I::GxBIterator{<:Any, <:Any, false, NeighborIterator()}) = increment(_rc_geti(I))
 @inline get_element(I::GxBIterator{<:Any, <:Any, true, NoIndexIterator()}) = getval(I)
 @inline get_element(::GxBIterator{<:Any, <:Any, false, NoIndexIterator()}) = throw(ArgumentError("Must iterate over either indices or values."))
+@inline get_element(I::GxBIterator{<:Any, <:Any, <:Any, IteratorIterator()}) = I
+
 
 struct VectorIterator{B, O, T, IterateValues, IterationType, I<:GxBIterator}
     iterator::I
@@ -228,10 +242,14 @@ end
 end
 
 get_element(I::VectorIterator) = get_element(I.iterator)
+get_element(I::VectorIterator{<:Any, <:Any, <:Any, <:Any, IteratorIterator()}) = I
 
 const RowIterator{B, T, IterateValues, IterationType} = VectorIterator{B, RowMajor(), T, IterateValues, IterationType}
 const ColIterator{B, T, IterateValues, IterationType} = VectorIterator{B, ColMajor(), T, IterateValues, IterationType}
 
+defaultiteration(::Integer) = NeighborIterator()
+defaultiteration(::AbstractVector) = IndicesIterator()
+
 RowIterator(A::AbstractGBArray, v, iteratevalues::Bool, indexiteration::IndexIteratorType = NeighborIterator()) =
     storageorder(A) === RowMajor() ? VectorIterator{iteratevalues, indexiteration}(A, v) : 
     throw(ArgumentError("A is not in RowMajor() order. Row iteration is only supported on RowMajor AbstractGBArrays. Try setstorageorder[!]"))
@@ -255,4 +273,9 @@ iteraterows(A::AbstractGBArray, v, indexiteration::IndexIteratorType = NeighborI
 function Base.iterate(I::VectorIterator)
     return _seek(I.iterator, I.v isa Int64 ? I.v : I.v.start) == LibGraphBLAS.GrB_NO_VALUE ? knext(I) : (get_element(I), nothing)
 end
-Base.iterate(I::VectorIterator, ::Nothing) = inext(I)
+Base.iterate(I::VectorIterator, ::Nothing) = inext(I)
+
+Base.getindex(A::AbstractArray, I::GxBIterator) = getindex(A, getrow(I), getcol(I))
+Base.getindex(A::AbstractArray, v::VectorIterator) = getindex(A, v.iterator)
+Base.setindex!(A::AbstractArray, x, I::GxBIterator) = setindex!(A, x, getrow(I), getcol(I))
+Base.setindex!(A::AbstractArray, x, v::VectorIterator) = setindex!(A, x, v.iterator)

src/linalg.jl

@@ -0,0 +1,4 @@
+# not quite as fast as it *could* be, we have to construct the sparse representation of
+# a diagonal. Which means ~triple the size of a basic Diagonal matrix.
+# SSGrB doesn't have an internal representation for a Diagonal.
+LinearAlgebra.:\(D::Diagonal, B::AbstractGBMatrix) = *(D, B, (any, \))

src/operations/broadcasts.jl

@@ -63,6 +63,11 @@ modifying(::typeof(*)) = mul!
 modifying(::typeof(eadd)) = eadd!
 modifying(::typeof(emul)) = emul!
 
+# TODO: Fix this horrifically ugly function.
+# TODO: Remove the requirement that vector -> matrix
+# broadcast be done on certain side.
+# We can add an `iscommutative` to solve this.
+# as well as some function to get the reversed version of the operator.
 @inline function Base.copy(bc::Broadcast.Broadcasted{GBMatrixStyle})
     f = bc.f
     l = length(bc.args)
@@ -92,6 +97,28 @@ modifying(::typeof(emul)) = emul!
         if right isa StridedArray
             right = pack(right; fill = left isa GBArrayOrTranspose ? getfill(left) : nothing)
         end
+        # TODO: We want to expose the broadcasting of Vectors into Matrices.
+        # The only problem is we need a notion of commutativity.
+        # This works fine for builtins, we can define commutativity pretty easily
+        # for many operations. But for non-builtins we'd need an API of sorts.
+        # To get around this for now we will require that Vectors be on the left
+        # and transposed vectors be on the right.
+        if left isa AbstractGBVector && right isa GBMatrixOrTranspose
+            return *(Diagonal(left), right, (any, f))
+        end
+        if left isa GBMatrixOrTranspose && right isa Transpose{<:Any, <:AbstractGBVector}
+            return *(left, Diagonal(right), (any, f))
+        end
+        if left isa GBMatrixOrTranspose && right isa AbstractGBVector
+            throw(ArgumentError(
+                "Broadcasting a GBVector into a GBMatrix is only currently " *
+                "supported with the GBVector on the left."))
+        end
+        if right isa GBMatrixOrTranspose && left isa Transpose{<:Any, <:AbstractGBVector}
+            throw(ArgumentError(
+                "Broadcasting a Transpose{<:Any, <:AbstractGBVector} into a GBMatrix" *
+                " is only currently supported with the GBVector on the right."))
+        end
         if left isa GBArrayOrTranspose && right isa GBArrayOrTranspose
             add = defaultadd(f)
             return add(left, right, f)

src/operations/extract.jl

@@ -68,10 +68,10 @@ function extract!(
     desc = _handledescriptor(desc; out=C, in1 = A)
     mask = _handlemask!(desc, mask)
     I = decrement!(I)
-    J = decrement!(J)
+    I !== J && (J = decrement!(J))
     @wraperror LibGraphBLAS.GrB_Matrix_extract(C, mask, _handleaccum(accum, storedeltype(C)), parent(A), I, ni, J, nj, desc)
-    I isa AbstractVector && increment!(I)
-    J isa AbstractVector && increment!(J)
+    increment!(I)
+    I !== J && increment!(J)
     return C
 end
 """
@@ -140,7 +140,7 @@ function extract!(
     desc = _handledescriptor(desc; out=w)
     mask = _handlemask!(desc, mask)
     @wraperror LibGraphBLAS.GrB_Matrix_extract(w, mask, _handleaccum(accum, storedeltype(w)), u, I, ni, UInt64[0], 1, desc)
-    I isa AbstractVector && increment!(I)
+    increment!(I)
     return w
 end
 

src/operations/mul.jl

@@ -191,13 +191,13 @@ function Base.:*(rig::TypedSemiring)
 end
 
 # Diagonal
-function Base.:*(D::Diagonal, A::G) where 
+function Base.:*(D::Diagonal, A::G, op = (+, *); mask = nothing, accum = nothing, desc = nothing) where 
     {G <: Union{Transpose{T, <:SuiteSparseGraphBLAS.AbstractGBArray{T1, F, O}} where {T, T1, F, O}, 
         SuiteSparseGraphBLAS.AbstractGBArray{T, F, O, 2} where {T, F, O}}}
-    return *(G(D), A)
+    return *(G(D), A, op; mask, accum, desc)
 end
-function Base.:*(A::G, D::Diagonal) where 
+function Base.:*(A::G, D::Diagonal, op = (+, *); mask = nothing, accum = nothing, desc = nothing) where 
     {G <: Union{Transpose{T, <:SuiteSparseGraphBLAS.AbstractGBArray{T1, F, O}} where {T, T1, F, O}, 
         SuiteSparseGraphBLAS.AbstractGBArray{T, F, O, 2} where {T, F, O}}}
-    return *(G(D), A)
+    return *(A, G(D), op; mask, accum, desc)
 end