diff --git a/Project.toml b/Project.toml index 709307bb..503b950e 100644 --- a/Project.toml +++ b/Project.toml @@ -5,7 +5,6 @@ version = "0.7.2" [deps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" -HyperSparseMatrices = "c7efdb1c-7caa-4c7d-9b5e-9093f9323c7c" KLU = "ef3ab10e-7fda-4108-b977-705223b18434" Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -21,7 +20,6 @@ SuiteSparse = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9" [compat] ChainRulesCore = "1" -HyperSparseMatrices = "0.2" MacroTools = "0.5" Preferences = "1" SSGraphBLAS_jll = "6.2.1" diff --git a/docs/src/binaryops.md b/docs/src/binaryops.md index 8498a504..3600ba16 100644 --- a/docs/src/binaryops.md +++ b/docs/src/binaryops.md @@ -27,9 +27,7 @@ Internally functions are lowered like this: ```@repl using SuiteSparseGraphBLAS -op = BinaryOp(+) - -typedop = op(Int64, Int64) +typedop = binaryop(+, Int64, Int64) eadd(GBVector([1,2]), GBVector([3,4]), typedop) ``` diff --git a/docs/src/operators.md b/docs/src/operators.md index d4b88660..48e4c0df 100644 --- a/docs/src/operators.md +++ b/docs/src/operators.md @@ -51,8 +51,7 @@ Operators are lowered from a Julia function to a container like `BinaryOp` or `S using SuiteSparseGraphBLAS ``` ```@repl operators -b = BinaryOp(+) -b(Int32) +b = binaryop(+, Int32) s = Semiring(max, +) s(Float64) diff --git a/docs/src/udfs.md b/docs/src/udfs.md index 3b355006..a6ac173d 100644 --- a/docs/src/udfs.md +++ b/docs/src/udfs.md @@ -7,8 +7,8 @@ GraphBLAS supports users to supply functions as operators. Constructors exported - `UnaryOp(name::String, fn::Function, [type | types | ztype, xtype | ztypes, xtypes])` - `BinaryOp(name::String, fn::Function, [type | types | ztype, xtype | ztypes, xtypes])` -- `Monoid(name::String, binop::Union{AbstractBinaryOp, GrB_BinaryOp}, id::T, terminal::T = nothing)`: all types must be the same. -- `Semiring(name::String, add::[GrB_Monoid | AbstractMonoid], mul::[GrB_BinaryOp | AbstractBinaryOp])` +- `Monoid(name::String, binop::Union{GrB_BinaryOp}, id::T, terminal::T = nothing)`: all types must be the same. +- `Semiring(name::String, add::[GrB_Monoid | AbstractMonoid], mul::GrB_BinaryOp)` `GrB_` prefixed arguments are typed operators, such as the result of `UnaryOps.COS[Float64]`. Type arguments may be single types or vectors of types. diff --git a/docs/src/unaryops.md b/docs/src/unaryops.md index 3b63bd28..85f1333e 100644 --- a/docs/src/unaryops.md +++ b/docs/src/unaryops.md @@ -19,9 +19,7 @@ Internally functions are lowered like this: ```@repl using SuiteSparseGraphBLAS -op = UnaryOp(sin) - -typedop = op(Float64) +op = unaryop(sin, Float64) map(typedop, GBVector([1.5, 0, pi])) ``` diff --git a/src/SuiteSparseGraphBLAS.jl b/src/SuiteSparseGraphBLAS.jl index ba02ecc3..00e5ddb7 100644 --- a/src/SuiteSparseGraphBLAS.jl +++ b/src/SuiteSparseGraphBLAS.jl @@ -25,8 +25,6 @@ using Serialization using StorageOrders export ColMajor, RowMajor, storageorder #reexports from StorageOrders - -using HyperSparseMatrices include("abstracts.jl") include("libutils.jl") @@ -101,7 +99,7 @@ include("oriented.jl") export SparseArrayCompat export LibGraphBLAS # export UnaryOps, BinaryOps, Monoids, Semirings #Submodules -export UnaryOp, BinaryOp, Monoid, Semiring #UDFs +export unaryop, binaryop, Monoid, semiring #UDFs export Descriptor #Types export gbset, gbget # global and object specific options. # export xtype, ytype, ztype #Determine input/output types of operators diff --git a/src/abstractgbarray.jl b/src/abstractgbarray.jl index f62fdaa6..2e2559bc 100644 --- a/src/abstractgbarray.jl +++ b/src/abstractgbarray.jl @@ -1,11 +1,11 @@ # AbstractGBArray functions: -function SparseArrays.nnz(A::AbsGBArrayOrTranspose) +function SparseArrays.nnz(A::GBArrayOrTranspose) nvals = Ref{LibGraphBLAS.GrB_Index}() @wraperror LibGraphBLAS.GrB_Matrix_nvals(nvals, gbpointer(parent(A))) return Int64(nvals[]) end -Base.eltype(::Type{AbstractGBArray{T}}) where{T} = T +Base.eltype(::Type{GBArrayOrTranspose{T}}) where{T} = T """ empty!(v::GBVector) @@ -14,32 +14,32 @@ Base.eltype(::Type{AbstractGBArray{T}}) where{T} = T Clear all the entries from the GBArray. Does not modify the type or dimensions. """ -function Base.empty!(A::AbsGBArrayOrTranspose) +function Base.empty!(A::GBArrayOrTranspose) @wraperror LibGraphBLAS.GrB_Matrix_clear(gbpointer(parent(A))) return A end -function Base.Matrix(A::AbstractGBMatrix) +function Base.Matrix(A::GBArrayOrTranspose) sparsity = sparsitystatus(A) T = copy(A) # We copy here to 1. avoid densifying A, and 2. to avoid destroying A. return unpack!(T, Dense()) end -function Base.Vector(v::AbstractGBVector) +function Base.Vector(v::GBVectorOrTranspose) sparsity = sparsitystatus(v) T = copy(v) # avoid densifying v and destroying v. return unpack!(T, Dense()) end -function SparseArrays.SparseMatrixCSC(A::AbstractGBArray) +function SparseArrays.SparseMatrixCSC(A::GBArrayOrTranspose) sparsity = sparsitystatus(A) T = copy(A) # avoid changing sparsity of A and destroying it. return unpack!(T, SparseMatrixCSC) end -function SparseArrays.SparseVector(v::AbstractGBVector) +function SparseArrays.SparseVector(v::GBVectorOrTranspose) sparsity = sparsitystatus(v) - T = copy(A) # avoid changing sparsity of v and destroying it. + T = copy(v) # avoid changing sparsity of v and destroying it. return unpack!(T, SparseVector) end @@ -94,7 +94,7 @@ for T ∈ valid_vec function build(A::AbstractGBMatrix{$T}, I::AbstractVector{<:Integer}, J::AbstractVector{<:Integer}, X::AbstractVector{$T}; combine = + ) - combine = BinaryOp(combine)($T) + combine = binaryop(combine, $T) I isa Vector || (I = collect(I)) J isa Vector || (J = collect(J)) X isa Vector || (X = collect(X)) @@ -181,7 +181,7 @@ function build( A::AbstractGBMatrix{T}, I::AbstractVector{<:Integer}, J::AbstractVector{<:Integer}, X::AbstractVector{T}; combine = + ) where {T} - combine = BinaryOp(combine)(T) + combine = binaryop(combine, T) I isa Vector || (I = collect(I)) J isa Vector || (J = collect(J)) X isa Vector || (X = collect(X)) @@ -314,7 +314,7 @@ Assign a submatrix of `A` to `C`. Equivalent to [`assign!`](@ref) except that # Keywords - `mask::Union{Nothing, GBMatrix} = nothing`: mask where `size(M) == size(A)`. -- `accum::Union{Nothing, Function, AbstractBinaryOp} = nothing`: binary accumulator operation +- `accum::Union{Nothing, Function} = nothing`: binary accumulator operation where `C[i,j] = accum(C[i,j], T[i,j])` where T is the result of this function before accum is applied. - `desc::Union{Nothing, Descriptor} = nothing` @@ -325,7 +325,7 @@ Assign a submatrix of `A` to `C`. Equivalent to [`assign!`](@ref) except that - `GrB_DIMENSION_MISMATCH`: If `size(A) != (max(I), max(J))` or `size(A) != size(mask)`. """ function subassign!( - C::AbstractGBArray, A::AbstractGBArray, I, J; + C::AbstractGBArray, A::GBArrayOrTranspose, I, J; mask = nothing, accum = nothing, desc = nothing ) I, ni = idx(I) @@ -397,7 +397,7 @@ Assign a submatrix of `A` to `C`. Equivalent to [`subassign!`](@ref) except that # Keywords - `mask::Union{Nothing, GBMatrix} = nothing`: mask where `size(M) == size(C)`. -- `accum::Union{Nothing, Function, AbstractBinaryOp} = nothing`: binary accumulator operation +- `accum::Union{Nothing, Function} = nothing`: binary accumulator operation where `C[i,j] = accum(C[i,j], T[i,j])` where T is the result of this function before accum is applied. - `desc::Union{Nothing, Descriptor} = nothing` @@ -408,7 +408,7 @@ Assign a submatrix of `A` to `C`. Equivalent to [`subassign!`](@ref) except that - `GrB_DIMENSION_MISMATCH`: If `size(A) != (max(I), max(J))` or `size(C) != size(mask)`. """ function assign!( - C::AbstractGBMatrix, A::AbstractGBVector, I, J; + C::AbstractGBMatrix, A::GBArrayOrTranspose, I, J; mask = nothing, accum = nothing, desc = nothing ) I, ni = idx(I) @@ -417,16 +417,16 @@ function assign!( I = decrement!(I) J = decrement!(J) # we know A isn't adjoint/transpose on input - desc = _handledescriptor(desc) - @wraperror LibGraphBLAS.GrB_Matrix_assign(gbpointer(C), mask, getaccum(accum, eltype(C)), gbpointer(A), I, ni, J, nj, desc) + desc = _handledescriptor(desc; in1=A) + @wraperror LibGraphBLAS.GrB_Matrix_assign(gbpointer(C), mask, getaccum(accum, eltype(C)), gbpointer(parent(A)), I, ni, J, nj, desc) increment!(I) increment!(J) return A end -function assign!(C::AbstractGBArray, x, I, J; +function assign!(C::AbstractGBArray{T}, x, I, J; mask = nothing, accum = nothing, desc = nothing -) +) where T x = typeof(x) === T ? x : convert(T, x) I, ni = idx(I) J, nj = idx(J) @@ -467,7 +467,7 @@ end Base.eltype(::Type{AbstractGBVector{T}}) where{T} = T function Base.deleteat!(v::AbstractGBVector, i) - @wraperror LibGraphBLAS.GrB_Matrix_removeElement(gbpointer(v), decrement!(i), 1) + @wraperror LibGraphBLAS.GrB_Matrix_removeElement(gbpointer(v), decrement!(i), 0) return v end @@ -520,7 +520,7 @@ for T ∈ valid_vec I isa Vector || (I = collect(I)) X isa Vector || (X = collect(X)) length(X) == length(I) || DimensionMismatch("I and X must have the same length") - combine = BinaryOp(combine)($T) + combine = binaryop(combine, $T) decrement!(I) @wraperror LibGraphBLAS.$func( Ptr{LibGraphBLAS.GrB_Vector}(gbpointer(v)), @@ -606,7 +606,7 @@ function build(v::AbstractGBVector{T}, I::Vector{<:Integer}, X::Vector{T}; combi I isa Vector || (I = collect(I)) X isa Vector || (X = collect(X)) length(X) == length(I) || DimensionMismatch("I and X must have the same length") - combine = BinaryOp(combine)(T) + combine = binaryop(combine, T) decrement!(I) @wraperror LibGraphBLAS.GrB_Matrix_build_UDT( Ptr{LibGraphBLAS.GrB_Vector}(gbpointer(v)), diff --git a/src/abstracts.jl b/src/abstracts.jl index f9d38a35..470b8ab4 100644 --- a/src/abstracts.jl +++ b/src/abstracts.jl @@ -1,11 +1,8 @@ abstract type AbstractGBType end abstract type AbstractDescriptor end abstract type AbstractOp end -abstract type AbstractUnaryOp <: AbstractOp end -abstract type AbstractBinaryOp <: AbstractOp end abstract type AbstractSelectOp <: AbstractOp end abstract type AbstractMonoid <: AbstractOp end -abstract type AbstractSemiring <: AbstractOp end abstract type AbstractTypedOp{Z} end abstract type AbstractGBArray{T, N, F} <: AbstractSparseArray{T, UInt64, N} end diff --git a/src/chainrules/chainruleutils.jl b/src/chainrules/chainruleutils.jl index cb71c792..10d19a1c 100644 --- a/src/chainrules/chainruleutils.jl +++ b/src/chainrules/chainruleutils.jl @@ -4,4 +4,4 @@ using ChainRulesCore const RealOrComplex = Union{Real, Complex} # LinearAlgebra.norm doesn't like the nothings. -LinearAlgebra.norm(A::GBArray, p::Real=2) = norm(nonzeros(A), p) +LinearAlgebra.norm(A::GBVecOrMat, p::Real=2) = norm(nonzeros(A), p) diff --git a/src/chainrules/ewiserules.jl b/src/chainrules/ewiserules.jl index 562bbc52..bdddae21 100644 --- a/src/chainrules/ewiserules.jl +++ b/src/chainrules/ewiserules.jl @@ -2,19 +2,19 @@ function frule( (_, ΔA, ΔB, _)::Tuple, ::typeof(emul), - A::GBArray, - B::GBArray, + A::AbstractGBArray, + B::AbstractGBArray, ::typeof(*) ) Ω = emul(A, B, *) ∂Ω = emul(unthunk(ΔA), B, *) + emul(unthunk(ΔB), A, *) return Ω, ∂Ω end -function frule((_, ΔA, ΔB)::Tuple, ::typeof(emul), A::GBArray, B::GBArray) +function frule((_, ΔA, ΔB)::Tuple, ::typeof(emul), A::AbstractGBArray, B::AbstractGBArray) return frule((nothing, ΔA, ΔB, nothing), emul, A, B, *) end -function rrule(::typeof(emul), A::GBArray, B::GBArray, ::typeof(*)) +function rrule(::typeof(emul), A::AbstractGBArray, B::AbstractGBArray, ::typeof(*)) function timespullback(ΔΩ) ∂A = emul(unthunk(ΔΩ), B) ∂B = emul(unthunk(ΔΩ), A) @@ -23,7 +23,7 @@ function rrule(::typeof(emul), A::GBArray, B::GBArray, ::typeof(*)) return emul(A, B, *), timespullback end -function rrule(::typeof(emul), A::GBArray, B::GBArray) +function rrule(::typeof(emul), A::AbstractGBArray, B::AbstractGBArray) Ω, fullpb = rrule(emul, A, B, *) emulpb(ΔΩ) = fullpb(ΔΩ)[1:3] return Ω, emulpb @@ -39,19 +39,19 @@ end function frule( (_, ΔA, ΔB, _)::Tuple, ::typeof(eadd), - A::GBArray, - B::GBArray, + A::AbstractGBArray, + B::AbstractGBArray, ::typeof(+) ) Ω = eadd(A, B, +) ∂Ω = eadd(unthunk(ΔA), unthunk(ΔB), +) return Ω, ∂Ω end -function frule((_, ΔA, ΔB)::Tuple, ::typeof(eadd), A::GBArray, B::GBArray) +function frule((_, ΔA, ΔB)::Tuple, ::typeof(eadd), A::AbstractGBArray, B::AbstractGBArray) return frule((nothing, ΔA, ΔB, nothing), eadd, A, B, +) end -function rrule(::typeof(eadd), A::GBArray, B::GBArray, ::typeof(+)) +function rrule(::typeof(eadd), A::AbstractGBArray, B::AbstractGBArray, ::typeof(+)) function pluspullback(ΔΩ) return ( NoTangent(), @@ -63,7 +63,7 @@ function rrule(::typeof(eadd), A::GBArray, B::GBArray, ::typeof(+)) return eadd(A, B, +), pluspullback end -function rrule(::typeof(eadd), A::GBArray, B::GBArray) +function rrule(::typeof(eadd), A::AbstractGBArray, B::AbstractGBArray) Ω, fullpb = rrule(eadd, A, B, +) eaddpb(ΔΩ) = fullpb(ΔΩ)[1:3] return Ω, eaddpb diff --git a/src/chainrules/maprules.jl b/src/chainrules/maprules.jl index ef728aea..96de9cc9 100644 --- a/src/chainrules/maprules.jl +++ b/src/chainrules/maprules.jl @@ -21,7 +21,7 @@ macro scalarapplyrule(func, derivative) (_, _, $(esc(:ΔA)))::Tuple, ::typeof(apply), ::typeof($(func)), - $(esc(:A))::GBArray + $(esc(:A))::AbstractGBArray ) $(esc(:Ω)) = apply($(esc(func)), $(esc(:A))) return $(esc(:Ω)), $(esc(derivative)) .* unthunk($(esc(:ΔA))) @@ -29,7 +29,7 @@ macro scalarapplyrule(func, derivative) function ChainRulesCore.rrule( ::typeof(apply), ::typeof($(func)), - $(esc(:A))::GBArray + $(esc(:A))::AbstractGBArray ) $(esc(:Ω)) = apply($(esc(func)), $(esc(:A))) function applyback($(esc(:ΔA))) @@ -75,10 +75,10 @@ function frule( (_, _, ΔA)::Tuple, ::typeof(apply), ::typeof(identity), - A::GBArray + A::AbstractGBArray ) return (A, ΔA) end -function rrule(::typeof(apply), ::typeof(identity), A::GBArray) +function rrule(::typeof(apply), ::typeof(identity), A::AbstractGBArray) return A, (ΔΩ) -> (NoTangent(), NoTangent(), ΔΩ) end diff --git a/src/chainrules/mulrules.jl b/src/chainrules/mulrules.jl index 766f68d0..81d2e404 100644 --- a/src/chainrules/mulrules.jl +++ b/src/chainrules/mulrules.jl @@ -3,16 +3,16 @@ function frule( (_, ΔA, ΔB)::Tuple, ::typeof(*), - A::GBArray, - B::GBArray + A::AbstractGBArray, + B::AbstractGBArray ) frule((nothing, ΔA, ΔB, nothing), *, A, B, (+, *)) end function frule( (_, ΔA, ΔB, _)::Tuple, ::typeof(*), - A::GBArray, - B::GBArray, + A::AbstractGBArray, + B::AbstractGBArray, ::typeof((+, *)) ) Ω = *(A, B, (+, *)) @@ -22,8 +22,8 @@ end function rrule( ::typeof(*), - A::GBArray, - B::GBArray, + A::AbstractGBArray, + B::AbstractGBArray, ::typeof((+, *)) ) function mulpullback(ΔΩ) @@ -37,8 +37,8 @@ end function rrule( ::typeof(*), - A::GBArray, - B::GBArray + A::AbstractGBArray, + B::AbstractGBArray ) Ω, mulpullback = rrule(*, A, B, (+, *)) pullback(ΔΩ) = mulpullback(ΔΩ)[1:3] @@ -50,8 +50,8 @@ end # Missing frule here. function rrule( ::typeof(*), - A::GBArray, - B::GBArray, + A::AbstractGBArray, + B::AbstractGBArray, ::typeof((+, /)) ) function mulpullback(ΔΩ) @@ -66,8 +66,8 @@ end function frule( (_, ΔA, ΔB, _)::Tuple, ::typeof(*), - A::GBArray, - B::GBArray, + A::AbstractGBArray, + B::AbstractGBArray, ::typeof((+, +)) ) Ω = *(A, B, (+, +)) @@ -77,8 +77,8 @@ end function rrule( ::typeof(*), - A::GBArray, - B::GBArray, + A::AbstractGBArray, + B::AbstractGBArray, ::typeof((+, +)) ) function mulpullback(ΔΩ) @@ -93,8 +93,8 @@ end function frule( (_, ΔA, ΔB, _)::Tuple, ::typeof(*), - A::GBArray, - B::GBArray, + A::AbstractGBArray, + B::AbstractGBArray, ::typeof((+, -)) ) Ω = *(A, B, (+, -)) @@ -104,8 +104,8 @@ end function rrule( ::typeof(*), - A::GBArray, - B::GBArray, + A::AbstractGBArray, + B::AbstractGBArray, ::typeof((+, -)) ) function mulpullback(ΔΩ) @@ -120,8 +120,8 @@ end function frule( (_, ΔA, ΔB, _)::Tuple, ::typeof(*), - A::GBArray, - B::GBArray, + A::AbstractGBArray, + B::AbstractGBArray, ::typeof((+, first)) ) Ω = *(A, B, (+, first)) @@ -131,8 +131,8 @@ end function rrule( ::typeof(*), - A::GBArray, - B::GBArray, + A::AbstractGBArray, + B::AbstractGBArray, ::typeof((+, first)) ) function mulpullback(ΔΩ) @@ -147,8 +147,8 @@ end function frule( (_, ΔA, ΔB, _)::Tuple, ::typeof(*), - A::GBArray, - B::GBArray, + A::AbstractGBArray, + B::AbstractGBArray, ::typeof((+, second)) ) Ω = *(A, B, (+, second)) @@ -158,8 +158,8 @@ end function rrule( ::typeof(*), - A::GBArray, - B::GBArray, + A::AbstractGBArray, + B::AbstractGBArray, ::typeof((+, second)) ) function mulpullback(ΔΩ) diff --git a/src/chainrules/selectrules.jl b/src/chainrules/selectrules.jl index d3d2d198..db2a7f25 100644 --- a/src/chainrules/selectrules.jl +++ b/src/chainrules/selectrules.jl @@ -3,7 +3,7 @@ function frule( (_, _, ΔA)::Tuple, ::typeof(select), op::Union{Function, SelectUnion}, - A::GBArray + A::AbstractGBArray ) Ω = select(op, A) ∂Ω = mask(unthunk(ΔA), Ω, structural = true) @@ -15,7 +15,7 @@ function frule( (_, _, ΔA, _)::Tuple, ::typeof(select), op::Union{Function, SelectUnion}, - A::GBArray, + A::AbstractGBArray, thunk::Union{GBScalar, Nothing, valid_union} ) Ω = select(op, A, thunk) @@ -26,7 +26,7 @@ end function rrule( ::typeof(select), op::Union{Function, SelectUnion}, - A::GBArray + A::AbstractGBArray ) out = select(op, A) function selectback(ΔΩ) @@ -39,7 +39,7 @@ end function rrule( ::typeof(select), op::Union{Function, SelectUnion}, - A::GBArray, + A::AbstractGBArray, thunk::Union{GBScalar, Nothing, valid_union} ) out = select(op, A, thunk) diff --git a/src/constants.jl b/src/constants.jl index e8318ecb..7448ad14 100644 --- a/src/constants.jl +++ b/src/constants.jl @@ -1,11 +1,7 @@ const GBVecOrMat{T} = Union{<:AbstractGBVector{T}, <:AbstractGBMatrix{T}} -const GBMatOrTranspose{T} = Union{<:AbstractGBMatrix{T}, Transpose{T, <:AbstractGBMatrix{T}}} -const GBVecOrTranspose{T} = Union{<:AbstractGBVector{T}, Transpose{T, <:AbstractGBVector{T}}} -const GBArray{T} = Union{<:GBVecOrTranspose{T}, <:GBMatOrTranspose{T}} - -const GBMatrixOrTranspose{T} = Union{<:GBMatrix{T}, Transpose{T, <:GBMatrix{T}}} -const GBVectorOrTranspose{T} = Union{<:GBVector{T}, Transpose{T, <:GBVector{T}}} -const AbsGBArrayOrTranspose{T} = Union{<:AbstractGBArray{T}, Transpose{T, <:AbstractGBArray{T}}} +const GBMatrixOrTranspose{T} = Union{<:AbstractGBMatrix{T}, Transpose{T, <:AbstractGBMatrix{T}}} +const GBVectorOrTranspose{T} = Union{<:AbstractGBVector{T}, Transpose{T, <:AbstractGBVector{T}}} +const GBArrayOrTranspose{T} = Union{<:AbstractGBArray{T}, Transpose{T, <:AbstractGBArray{T}}} const ptrtogbtype = IdDict{Ptr, GBType}() @@ -28,8 +24,6 @@ const MonoidBinaryOrRig = Union{ TypedMonoid, TypedSemiring, TypedBinaryOperator, - AbstractSemiring, - AbstractBinaryOp, AbstractMonoid } diff --git a/src/gbtypes.jl b/src/gbtypes.jl index 75d70c95..bf4d13ef 100644 --- a/src/gbtypes.jl +++ b/src/gbtypes.jl @@ -47,6 +47,10 @@ end gbpointer(T::GBType) = T.p + + +const GBTYPES = IdDict{DataType, GBType}() + """ gbtype(x) @@ -54,9 +58,14 @@ Determine the GBType equivalent of a Julia primitive type. See also: [`juliatype`](@ref) """ -function gbtype end +function gbtype(::Type{T}; builtin = false, loaded = false, typestr = string(T)) where T + (get!(GBTYPES, T) do + return GBType{T}(builtin, loaded, LibGraphBLAS.GrB_Type(), typestr) + end)::GBType{T} +end -macro gbtype(expr...) +# For maybe a tiny bit of additional speed we'll overload `gbtype` for builtins. +macro _gbtype(expr...) jtype = expr[1] if length(expr) == 2 @@ -112,17 +121,17 @@ See also: [`gbtype`](@ref) """ juliatype(::GBType{T}) where {T} = T -@gbtype Bool GrB_BOOL -@gbtype Int8 GrB_INT8 -@gbtype UInt8 GrB_UINT8 -@gbtype Int16 GrB_INT16 -@gbtype UInt16 GrB_UINT16 -@gbtype Int32 GrB_INT32 -@gbtype UInt32 GrB_UINT32 -@gbtype Int64 GrB_INT64 -@gbtype UInt64 GrB_UINT64 -@gbtype Float32 GrB_FP32 -@gbtype Float64 GrB_FP64 -@gbtype ComplexF32 GxB_FC32 -@gbtype ComplexF64 GxB_FC64 +@_gbtype Bool GrB_BOOL +@_gbtype Int8 GrB_INT8 +@_gbtype UInt8 GrB_UINT8 +@_gbtype Int16 GrB_INT16 +@_gbtype UInt16 GrB_UINT16 +@_gbtype Int32 GrB_INT32 +@_gbtype UInt32 GrB_UINT32 +@_gbtype Int64 GrB_INT64 +@_gbtype UInt64 GrB_UINT64 +@_gbtype Float32 GrB_FP32 +@_gbtype Float64 GrB_FP64 +@_gbtype ComplexF32 GxB_FC32 +@_gbtype ComplexF64 GxB_FC64 diff --git a/src/matrix.jl b/src/matrix.jl index 0900d16f..67be931d 100644 --- a/src/matrix.jl +++ b/src/matrix.jl @@ -159,6 +159,6 @@ end #################### # Linear indexing -function Base.getindex(A::GBMatOrTranspose, v::AbstractVector) +function Base.getindex(A::GBMatrixOrTranspose, v::AbstractVector) throw("Not implemented") end diff --git a/src/operations/broadcasts.jl b/src/operations/broadcasts.jl index 7a70875f..95827cb0 100644 --- a/src/operations/broadcasts.jl +++ b/src/operations/broadcasts.jl @@ -80,7 +80,7 @@ modifying(::typeof(emul)) = emul! if right isa Broadcast.Broadcasted right = copy(right) end - if left isa GBArray && right isa GBArray + if left isa GBArrayOrTranspose && right isa GBArrayOrTranspose add = defaultadd(f) return add(left, right, f) else @@ -91,7 +91,7 @@ end mutatingop(::typeof(emul)) = emul! mutatingop(::typeof(eadd)) = eadd! mutatingop(::typeof(apply)) = apply! -@inline function Base.copyto!(C::AbsGBArrayOrTranspose, bc::Broadcast.Broadcasted{GBMatrixStyle}) +@inline function Base.copyto!(C::GBArrayOrTranspose, bc::Broadcast.Broadcasted{GBMatrixStyle}) l = length(bc.args) if l == 1 x = first(bc.args) @@ -139,7 +139,7 @@ mutatingop(::typeof(apply)) = apply! # If they're further nested broadcasts we can't fuse them, so just copy. subargleft isa Broadcast.Broadcasted && (subargleft = copy(subargleft)) subargright isa Broadcast.Broadcasted && (subargright = copy(subargright)) - if subargleft isa GBArray && subargright isa GBArray + if subargleft isa GBArrayOrTranspose && subargright isa GBArrayOrTranspose add = mutatingop(defaultadd(f)) return add(C, subargleft, subargright, f; accum) else @@ -156,7 +156,7 @@ mutatingop(::typeof(apply)) = apply! if right isa Broadcast.Broadcasted right = copy(right) end - if left isa GBArray && right isa GBArray + if left isa GBArrayOrTranspose && right isa GBArrayOrTranspose add = mutatingop(defaultadd(f)) return add(C, left, right, f) else @@ -189,7 +189,7 @@ end if right isa Broadcast.Broadcasted right = copy(right) end - if left isa GBArray && right isa GBArray + if left isa GBArrayOrTranspose && right isa GBArrayOrTranspose add = defaultadd(f) return add(left, right, f) else @@ -246,7 +246,7 @@ end # If they're further nested broadcasts we can't fuse them, so just copy. subargleft isa Broadcast.Broadcasted && (subargleft = copy(subargleft)) subargright isa Broadcast.Broadcasted && (subargright = copy(subargright)) - if subargleft isa GBArray && subargright isa GBArray + if subargleft isa GBArrayOrTranspose && subargright isa GBArrayOrTranspose add = mutatingop(defaultadd(f)) return add(C, subargleft, subargright, f; accum) else @@ -263,7 +263,7 @@ end if right isa Broadcast.Broadcasted right = copy(right) end - if left isa GBArray && right isa GBArray + if left isa GBArrayOrTranspose && right isa GBArrayOrTranspose add = mutatingop(defaultadd(f)) return add(C, left, right, f) else diff --git a/src/operations/concat.jl b/src/operations/concat.jl index f7f35597..1b625dcc 100644 --- a/src/operations/concat.jl +++ b/src/operations/concat.jl @@ -1,4 +1,4 @@ -function cat!(C::GBArray, tiles::AbstractArray{T}) where {T<:GBArray} +function cat!(C::GBVecOrMat, tiles::AbstractArray{T}) where {T<:GBVecOrMat} tiles = permutedims(tiles) @wraperror LibGraphBLAS.GxB_Matrix_concat(gbpointer(C), gbpointer.(tiles), size(tiles,2), size(tiles,1), C_NULL) return C @@ -32,10 +32,10 @@ function Base.cat(tiles::VecOrMat{<:Union{AbstractGBMatrix{T, F}, AbstractGBVect return cat!(C, tiles) end -vcat!(C, A::GBArray...) = cat!(C, collect(A)) -Base.vcat(A::GBArray...) = cat(collect(A)) +vcat!(C, A::GBArrayOrTranspose...) = cat!(C, collect(A)) +Base.vcat(A::GBArrayOrTranspose...) = cat(collect(A)) -hcat!(C, A::GBArray...) = cat!(C, permutedims(collect(A))) -Base.hcat(A::GBArray...) = cat(permutedims(collect(A))) +hcat!(C, A::GBArrayOrTranspose...) = cat!(C, permutedims(collect(A))) +Base.hcat(A::GBArrayOrTranspose...) = cat(permutedims(collect(A))) # TODO split. I don't necessarily see a great need for split though. We have indexing/slicing. \ No newline at end of file diff --git a/src/operations/ewise.jl b/src/operations/ewise.jl index 2dc38b4e..25eccc5b 100644 --- a/src/operations/ewise.jl +++ b/src/operations/ewise.jl @@ -1,5 +1,5 @@ """ - emul!(C::GBArray, A::GBArray, B::GBArray, op = *; kwargs...)::GBArray + emul!(C::GBArrayOrTranspose, A::GBArrayOrTranspose, B::GBArrayOrTranspose, op = *; kwargs...)::GBArrayOrTranspose Apply the binary operator `op` elementwise on the set intersection of `A` and `B`. Store or accumulate the result into C. When `op = *` this is equivalent to `A .* B`, @@ -9,22 +9,22 @@ The pattern of the result is the set intersection of `A` and `B`. For a set union equivalent see [`eadd!`](@ref). # Arguments -- `C::GBArray`: the output vector or matrix. -- `A, B::GBArray`: A GBVector or GBMatrix, possibly transposed. -- `op::Union{Function, AbstractBinaryOp, Monoid} = *`: the binary operation which is +- `C::GBArrayOrTranspose`: the output vector or matrix. +- `A, B::GBArrayOrTranspose`: A GBVector or GBMatrix, possibly transposed. +- `op::Union{Function, Monoid} = *`: the binary operation which is applied such that `C[i,j] = op(A[i,j], B[i,j])` for all `i,j` present in both `A` and `B`. # Keywords - `mask::Union{Nothing, GBVecOrMat} = nothing`: optional mask. -- `accum::Union{Nothing, Function, AbstractBinaryOp} = nothing`: binary accumulator operation +- `accum::Union{Nothing, Function} = nothing`: binary accumulator operation such that `C[i,j] = accum(C[i,j], T[i,j])` where T is the result of this function before accum is applied. - `desc::Union{Nothing, Descriptor} = nothing` """ function emul!( C::GBVecOrMat, - A::GBArray, - B::GBArray, + A::GBArrayOrTranspose, + B::GBArrayOrTranspose, op = *; mask = nothing, accum = nothing, @@ -33,7 +33,7 @@ function emul!( mask, accum = _handlenothings(mask, accum) desc = _handledescriptor(desc; in1=A, in2=B) size(C) == size(A) == size(B) || throw(DimensionMismatch()) - op = BinaryOp(op)(eltype(A), eltype(B)) + op = binaryop(op, eltype(A), eltype(B)) accum = getaccum(accum, eltype(C)) if op isa TypedBinaryOperator @wraperror LibGraphBLAS.GrB_Matrix_eWiseMult_BinaryOp(gbpointer(C), mask, accum, op, gbpointer(parent(A)), gbpointer(parent(B)), desc) @@ -44,7 +44,7 @@ function emul!( end """ - emul(A::GBArray, B::GBArray, op = *; kwargs...)::GBMatrix + emul(A::GBArrayOrTranspose, B::GBArrayOrTranspose, op = *; kwargs...)::GBMatrix Apply the binary operator `op` elementwise on the set intersection of `A` and `B`. When `op = *` this is equivalent to `A .* B`, however any binary operator may be substituted. @@ -53,13 +53,13 @@ The pattern of the result is the set intersection of `A` and `B`. For a set union equivalent see [`eadd`](@ref). # Arguments -- `A, B::GBArray`: A GBVector or GBMatrix, possibly transposed. -- `op::Union{Function, AbstractBinaryOp, Monoid} = *`: the binary operation which is +- `A, B::GBArrayOrTranspose`: A GBVector or GBMatrix, possibly transposed. +- `op::Union{Function, Monoid} = *`: the binary operation which is applied such that `C[i,j] = op(A[i,j], B[i,j])` for all `i,j` present in both `A` and `B`. # Keywords - `mask::Union{Nothing, GBVecOrMat} = nothing`: optional mask. -- `accum::Union{Nothing, AbstractBinaryOp} = nothing`: binary accumulator operation +- `accum::Union{Nothing} = nothing`: binary accumulator operation where `C[i,j] = accum(C[i,j], T[i,j])` where T is the result of this function before accum is applied. - `desc = nothing` @@ -68,8 +68,8 @@ union equivalent see [`eadd`](@ref). `A` and `B` or the binary operation if a type specific operation is provided. """ function emul( - A::GBArray, - B::GBArray, + A::GBArrayOrTranspose, + B::GBArrayOrTranspose, op = *; mask = nothing, accum = nothing, @@ -81,7 +81,7 @@ function emul( end """ - eadd!(C::GBVecOrMat, A::GBArray, B::GBArray, op = +; kwargs...)::GBVecOrMat + eadd!(C::GBVecOrMat, A::GBArrayOrTranspose, B::GBArrayOrTranspose, op = +; kwargs...)::GBVecOrMat Apply the binary operator `op` elementwise on the set union of `A` and `B`. Store or accumulate the result into C. When `op = +` this is equivalent to `A .+ B`, @@ -94,21 +94,21 @@ is an implicit zero returns `A[i,j]` **not** `A[i,j] op zero(T)`. For a set intersection equivalent see [`emul!`](@ref). # Arguments -- `C::GBArray`: the output vector or matrix. -- `A, B::GBArray`: A GBVector or GBMatrix, possibly transposed. -- `op::Union{Function, AbstractBinaryOp, Monoid} = +`: the binary operation which is +- `C::GBArrayOrTranspose`: the output vector or matrix. +- `A, B::GBArrayOrTranspose`: A GBVector or GBMatrix, possibly transposed. +- `op::Union{Function, Monoid} = +`: the binary operation which is applied such that `C[i,j] = op(A[i,j], B[i,j])` for all `i,j` present in either `A` and `B`. # Keywords - `mask::Union{Nothing, GBVecOrMat} = nothing`: optional mask. -- `accum::Union{Nothing, Function, AbstractBinaryOp} = nothing`: binary accumulator operation +- `accum::Union{Nothing, Function} = nothing`: binary accumulator operation such that `C[i,j] = accum(C[i,j], T[i,j])` where T is the result of this function before accum is applied. - `desc::Union{Nothing, Descriptor} = nothing` """ function eadd!( C::GBVecOrMat, - A::GBArray, - B::GBArray, + A::GBArrayOrTranspose, + B::GBArrayOrTranspose, op = +; mask = nothing, accum = nothing, @@ -117,18 +117,18 @@ function eadd!( mask, accum = _handlenothings(mask, accum) desc = _handledescriptor(desc; in1=A, in2 = B) size(C) == size(A) == size(B) || throw(DimensionMismatch()) - op = BinaryOp(op)(eltype(A), eltype(B)) + op = binaryop(op, eltype(A), eltype(B)) accum = getaccum(accum, eltype(C)) if op isa TypedBinaryOperator @wraperror LibGraphBLAS.GrB_Matrix_eWiseAdd_BinaryOp(gbpointer(C), mask, accum, op, gbpointer(parent(A)), gbpointer(parent(B)), desc) return C else - throw(ArgumentError("$op is not a valid monoid binary op or semiring.")) + throw(ArgumentError("$op is not a valid binary op.")) end end """ - eadd(A::GBArray, B::GBArray, op = +; kwargs...)::GBVecOrMat + eadd(A::GBArrayOrTranspose, B::GBArrayOrTranspose, op = +; kwargs...)::GBVecOrMat Apply the binary operator `op` elementwise on the set union of `A` and `B`. When `op = +` this is equivalent to `A .+ B`, however any binary operation may be substituted. @@ -140,19 +140,19 @@ is an implicit zero returns `A[i,j]` **not** `A[i,j] op zero(T)`. For a set intersection equivalent see [`emul`](@ref). # Arguments -- `A, B::GBArray`: A GBVector or GBMatrix, possibly transposed. -- `op::Union{Function, AbstractBinaryOp, Monoid} = +`: the binary operation which is +- `A, B::GBArrayOrTranspose`: A GBVector or GBMatrix, possibly transposed. +- `op::Union{Function, Monoid} = +`: the binary operation which is applied such that `C[i,j] = op(A[i,j], B[i,j])` for all `i,j` present in either `A` and `B`. # Keywords - `mask::Union{Nothing, GBVecOrMat} = nothing`: optional mask. -- `accum::Union{Nothing, Function, AbstractBinaryOp} = nothing`: binary accumulator operation +- `accum::Union{Nothing, Function} = nothing`: binary accumulator operation such that `C[i,j] = accum(C[i,j], T[i,j])` where T is the result of this function before accum is applied. - `desc::Union{Nothing, Descriptor} = nothing` """ function eadd( - A::GBArray, - B::GBArray, + A::GBArrayOrTranspose, + B::GBArrayOrTranspose, op = +; mask = nothing, accum = nothing, @@ -165,7 +165,7 @@ end """ - eunion!(C::GBVecOrMat, A::GBArray{T}, α::T B::GBArray, β::T, op = +; kwargs...)::GBVecOrMat + eunion!(C::GBVecOrMat, A::GBArrayOrTranspose{T}, α::T B::GBArrayOrTranspose, β::T, op = +; kwargs...)::GBVecOrMat Apply the binary operator `op` elementwise on the set union of `A` and `B`. Store or accumulate the result into C. When `op = +` this is equivalent to `A .+ B`, @@ -175,23 +175,23 @@ Unlike `eadd!` where an argument missing in `A` causes the `B` element to "pass- `eunion!` utilizes the `α` and `β` arguments for the missing operand elements. # Arguments -- `C::GBArray`: the output vector or matrix. -- `A, B::GBArray`: A GBVector or GBMatrix, possibly transposed. +- `C::GBArrayOrTranspose`: the output vector or matrix. +- `A, B::GBArrayOrTranspose`: A GBVector or GBMatrix, possibly transposed. - `α, β`: The fill-in value for `A` and `B` respectively. -- `op::Union{Function, AbstractBinaryOp, Monoid} = +`: the binary operation which is +- `op::Union{Function, Monoid} = +`: the binary operation which is applied such that `C[i,j] = op(A[i,j], B[i,j])` for all `i,j` present in either `A` and `B`. # Keywords - `mask::Union{Nothing, GBVecOrMat} = nothing`: optional mask. -- `accum::Union{Nothing, Function, AbstractBinaryOp} = nothing`: binary accumulator operation +- `accum::Union{Nothing, Function} = nothing`: binary accumulator operation such that `C[i,j] = accum(C[i,j], T[i,j])` where T is the result of this function before accum is applied. - `desc::Union{Nothing, Descriptor} = nothing` """ function eunion!( C::GBVecOrMat, - A::GBArray{T}, + A::GBArrayOrTranspose{T}, α::T, - B::GBArray{U}, + B::GBArrayOrTranspose{U}, β::U, op = +; mask = nothing, @@ -201,7 +201,7 @@ function eunion!( mask, accum = _handlenothings(mask, accum) desc = _handledescriptor(desc; in1=A, in2 = B) size(C) == size(A) == size(B) || throw(DimensionMismatch()) - op = BinaryOp(op)(eltype(A), eltype(B)) + op = binaryop(op, eltype(A), eltype(B)) accum = getaccum(accum, eltype(C)) if op isa TypedBinaryOperator @wraperror LibGraphBLAS.GxB_Matrix_eWiseUnion(gbpointer(C), mask, accum, op, gbpointer(parent(A)), GBScalar(α), gbpointer(parent(B)), GBScalar(β), desc) @@ -212,7 +212,7 @@ function eunion!( end """ - eunion(C::GBVecOrMat, A::GBArray{T}, α::T B::GBArray, β::T, op = +; kwargs...)::GBVecOrMat + eunion(C::GBVecOrMat, A::GBArrayOrTranspose{T}, α::T B::GBArrayOrTranspose, β::T, op = +; kwargs...)::GBVecOrMat Apply the binary operator `op` elementwise on the set union of `A` and `B`. When `op = +` this is equivalent to `A .+ B`, however any binary operation may be substituted. @@ -221,21 +221,21 @@ Unlike `eadd!` where an argument missing in `A` causes the `B` element to "pass- `eunion!` utilizes the `α` and `β` arguments for the missing operand elements. # Arguments -- `A, B::GBArray`: A GBVector or GBMatrix, possibly transposed. +- `A, B::GBArrayOrTranspose`: A GBVector or GBMatrix, possibly transposed. - `α, β`: The fill-in value for `A` and `B` respectively. -- `op::Union{Function, AbstractBinaryOp, Monoid} = +`: the binary operation which is +- `op::Union{Function, Monoid} = +`: the binary operation which is applied such that `C[i,j] = op(A[i,j], B[i,j])` for all `i,j` present in either `A` and `B`. # Keywords - `mask::Union{Nothing, GBVecOrMat} = nothing`: optional mask. -- `accum::Union{Nothing, Function, AbstractBinaryOp} = nothing`: binary accumulator operation +- `accum::Union{Nothing, Function} = nothing`: binary accumulator operation such that `C[i,j] = accum(C[i,j], T[i,j])` where T is the result of this function before accum is applied. - `desc::Union{Nothing, Descriptor} = nothing` """ function eunion( - A::GBArray{T}, + A::GBArrayOrTranspose{T}, α::T, - B::GBArray{U}, + B::GBArrayOrTranspose{U}, β::U, op = +; mask = nothing, @@ -247,11 +247,11 @@ function eunion( return eunion!(C, A, α, B, β, op; mask, accum, desc) end -function Base.:+(A::GBArray, B::GBArray) +function Base.:+(A::GBArrayOrTranspose, B::GBArrayOrTranspose) eadd(A, B, +) end -function Base.:-(A::GBArray, B::GBArray) +function Base.:-(A::GBArrayOrTranspose, B::GBArrayOrTranspose) eadd(A, B, -) end @@ -260,8 +260,8 @@ end ⊗(A, B, op; mask = nothing, accum = nothing, desc = nothing) = emul(A, B, op; mask, accum, desc) -⊕(f::Union{Function, BinaryUnion}) = (A, B; mask = nothing, accum = nothing, desc = nothing) -> +⊕(f::Union{Function, TypedBinaryOperator}) = (A, B; mask = nothing, accum = nothing, desc = nothing) -> eadd(A, B, f; mask, accum, desc) -⊗(f::Union{Function, BinaryUnion}) = (A, B; mask = nothing, accum = nothing, desc = nothing) -> +⊗(f::Union{Function, TypedBinaryOperator}) = (A, B; mask = nothing, accum = nothing, desc = nothing) -> emul(A, B, f; mask, accum, desc) diff --git a/src/operations/extract.jl b/src/operations/extract.jl index a2d2f334..e42281eb 100644 --- a/src/operations/extract.jl +++ b/src/operations/extract.jl @@ -30,7 +30,7 @@ function _outlength(u, I) end """ - extract!(C::GBMatrix, A::GBMatOrTranspose, I, J; kwargs...)::GBMatrix + extract!(C::GBMatrix, A::GBMatrixOrTranspose, I, J; kwargs...)::GBMatrix extract!(C::GBVector, A::GBVector, I; kwargs...)::GBVector Extract a submatrix or subvector from `A` into `C`. @@ -43,7 +43,7 @@ Extract a submatrix or subvector from `A` into `C`. # Keywords - `mask::Union{Nothing, GBArray} = nothing`: mask where `size(M) == (max(I), max(J))`. -- `accum::Union{Nothing, AbstractBinaryOp} = nothing`: binary accumulator operation +- `accum::Union{Nothing} = nothing`: binary accumulator operation where `C[i,j] = accum(C[i,j], T[i,j])` where T is the result of this function before accum is applied. - `desc::Union{Nothing, Descriptor} = nothing` @@ -57,7 +57,7 @@ Extract a submatrix or subvector from `A` into `C`. extract! function extract!( - C::AbstractGBMatrix, A::GBMatOrTranspose, I, J; + C::AbstractGBArray, A::GBMatrixOrTranspose, I, J; mask = nothing, accum = nothing, desc = nothing ) I, ni = idx(I) @@ -74,7 +74,7 @@ function extract!( return C end """ - extract(A::GBMatOrTranspose, I, J; kwargs...)::GBMatrix + extract(A::GBMatrixOrTranspose, I, J; kwargs...)::GBMatrix extract(A::GBVector, I; kwargs...)::GBVector Extract a submatrix or subvector from `A` @@ -86,7 +86,7 @@ Extract a submatrix or subvector from `A` # Keywords - `mask::Union{Nothing, GBArray} = nothing`: mask where `size(M) == (max(I), max(J))`. -- `accum::Union{Nothing, AbstractBinaryOp} = nothing`: binary accumulator operation +- `accum::Union{Nothing} = nothing`: binary accumulator operation where `C[i,j] = accum(C[i,j], T[i,j])` where T is the result of this function before accum is applied. - `desc::Descriptor = nothing` @@ -99,7 +99,7 @@ Extract a submatrix or subvector from `A` extract function extract( - A::GBMatOrTranspose, I, J; + A::GBMatrixOrTranspose, I, J; mask = nothing, accum = nothing, desc = nothing ) Ilen, Jlen = _outlength(A, I, J) @@ -107,6 +107,25 @@ function extract( return extract!(C, A, I, J; mask, accum, desc) end +function extract( + A::GBMatrixOrTranspose, ::Colon, J::Number; + mask = nothing, accum = nothing, desc = nothing +) + Ilen, _ = _outlength(A, :, J) + C = similar(A, Ilen) + return extract!(C, A, :, J; mask, accum, desc) +end + +# TODO: FINISH THIS WITH GxB_Matrix_reshape!!! +# function extract( +# A::GBMatrixOrTranspose, I::Number, ::Colon; +# mask = nothing, accum = nothing, desc = nothing +# ) +# _, Jlen = _outlength(A, I, :) +# C = similar(A, 1, Jlen) +# return extract!(C, A, I, :; mask, accum, desc) +# end + function extract!( w::AbstractGBVector, u::AbstractGBVector, I; mask = nothing, accum = nothing, desc = nothing @@ -115,7 +134,7 @@ function extract!( I = decrement!(I) desc = _handledescriptor(desc) mask === nothing && (mask = C_NULL) - @wraperror LibGraphBLAS.GrB_Matrix_extract(gbpointer(w), mask, getaccum(accum, eltype(w)), gbpointer(u), I, ni, UInt64[1], 1, desc) + @wraperror LibGraphBLAS.GrB_Matrix_extract(gbpointer(w), mask, getaccum(accum, eltype(w)), gbpointer(u), I, ni, UInt64[0], 1, desc) I isa Vector && increment!(I) return w end diff --git a/src/operations/kronecker.jl b/src/operations/kronecker.jl index b00b55c7..ae80d332 100644 --- a/src/operations/kronecker.jl +++ b/src/operations/kronecker.jl @@ -5,8 +5,8 @@ In-place version of [kron](@ref). """ function LinearAlgebra.kron!( C::GBVecOrMat, - A::GBArray, - B::GBArray, + A::GBArrayOrTranspose, + B::GBArrayOrTranspose, op = *; mask = nothing, accum = nothing, @@ -14,7 +14,7 @@ function LinearAlgebra.kron!( ) mask, accum = _handlenothings(mask, accum) desc = _handledescriptor(desc; in1=A, in2=B) - op = BinaryOp(op)(eltype(A), eltype(B)) + op = binaryop(op, eltype(A), eltype(B)) accum = getaccum(accum, eltype(C)) @wraperror LibGraphBLAS.GxB_kron(gbpointer(C), mask, accum, op, gbpointer(parent(A)), gbpointer(parent(B)), desc) return C @@ -34,13 +34,13 @@ Does not support `GBVector`s at this time. # Keywords - `mask::Union{Nothing, GBMatrix} = nothing`: optional mask. -- `accum::Union{Nothing, AbstractBinaryOp} = nothing`: binary accumulator operation +- `accum::Union{Nothing} = nothing`: binary accumulator operation where `C[i,j] = accum(C[i,j], T[i,j])` where T is the result of this function before accum is applied. - `desc = nothing` """ function LinearAlgebra.kron( - A::GBArray, - B::GBArray, + A::GBArrayOrTranspose, + B::GBArrayOrTranspose, op = *; mask = nothing, accum = nothing, diff --git a/src/operations/map.jl b/src/operations/map.jl index fec33101..42c265a6 100644 --- a/src/operations/map.jl +++ b/src/operations/map.jl @@ -1,26 +1,26 @@ function apply!( - op, C::GBVecOrMat, A::GBArray{T}; + op, C::GBVecOrMat, A::GBArrayOrTranspose{T}; mask = nothing, accum = nothing, desc = nothing ) where {T} mask, accum = _handlenothings(mask, accum) desc = _handledescriptor(desc; in1=A) - op = UnaryOp(op)(eltype(A)) + op = unaryop(op, eltype(A)) accum = getaccum(accum, eltype(C)) @wraperror LibGraphBLAS.GrB_Matrix_apply(gbpointer(C), mask, accum, op, gbpointer(parent(A)), desc) return C end function apply!( - op, A::GBArray{T}; + op, A::GBArrayOrTranspose{T}; mask = nothing, accum = nothing, desc = nothing ) where {T} return apply!(op, A, A; mask, accum, desc) end """ - apply(op::Union{Function, AbstractUnaryOp}, A::GBArray; kwargs...)::GBArray - apply(op::Union{Function, AbstractBinaryOp}, A::GBArray, x; kwargs...)::GBArray - apply(op::Union{Function, AbstractBinaryOp}, x, A::GBArray, kwargs...)::GBArray + apply(op::Union{Function, TypedUnaryOperator}, A::GBArrayOrTranspose; kwargs...)::GBArrayOrTranspose + apply(op::Union{Function}, A::GBArrayOrTranspose, x; kwargs...)::GBArrayOrTranspose + apply(op::Union{Function}, x, A::GBArrayOrTranspose, kwargs...)::GBArrayOrTranspose Transform a GBArray by applying `op` to each element. Equivalent to `Base.map` except for the additional `x` argument for mapping with a scalar. @@ -30,18 +30,18 @@ BinaryOps and two argument functions require the additional argument `x` which i substituted as the first or second operand of `op` depending on its position. # Arguments -- `op::Union{Function, AbstractUnaryOp, AbstractBinaryOp}` -- `A::GBArray` +- `op::Union{Function, TypedUnaryOperator}` +- `A::GBArrayOrTranspose` - `x`: Position dependent argument to binary operators. # Keywords - `mask::Union{Nothing, GBVecOrMat} = nothing`: optional mask. -- `accum::Union{Nothing, AbstractBinaryOp} = nothing`: binary accumulator operation +- `accum::Union{Nothing} = nothing`: binary accumulator operation where `C[i,j] = accum(C[i,j], T[i,j])` where T is the result of this function before accum is applied. - `desc::Union{Nothing, Descriptor} = nothing` """ function apply( - op, A::GBArray{T}; + op, A::GBArrayOrTranspose{T}; mask = nothing, accum = nothing, desc = nothing ) where {T} t = inferunarytype(eltype(A), op) @@ -49,26 +49,26 @@ function apply( end function apply!( - op, C::GBVecOrMat, x, A::GBArray{T}; + op, C::GBVecOrMat, x, A::GBArrayOrTranspose{T}; mask = nothing, accum = nothing, desc = nothing ) where {T} mask, accum = _handlenothings(mask, accum) desc = _handledescriptor(desc; in2=A) - op = BinaryOp(op)(eltype(A), typeof(x)) + op = binaryop(op, eltype(A), typeof(x)) accum = getaccum(accum, eltype(C)) @wraperror LibGraphBLAS.GxB_Matrix_apply_BinaryOp1st(gbpointer(C), mask, accum, op, GBScalar(x), gbpointer(parent(A)), desc) return C end function apply!( - op, x, A::GBArray{T}; + op, x, A::GBArrayOrTranspose{T}; mask = nothing, accum = nothing, desc = nothing ) where {T} return apply!(op, A, x, A; mask, accum, desc) end function apply( - op, x, A::GBArray{T}; + op, x, A::GBArrayOrTranspose{T}; mask = nothing, accum = nothing, desc = nothing ) where {T} t = inferbinarytype(typeof(x), eltype(A), op) @@ -76,55 +76,55 @@ function apply( end function apply!( - op, C::GBVecOrMat, A::GBArray{T}, x; + op, C::GBVecOrMat, A::GBArrayOrTranspose{T}, x; mask = nothing, accum = nothing, desc = nothing ) where {T} mask, accum = _handlenothings(mask, accum) desc = _handledescriptor(desc; in1=A) - op = BinaryOp(op)(eltype(A), typeof(x)) + op = binaryop(op, eltype(A), typeof(x)) accum = getaccum(accum, eltype(C)) @wraperror LibGraphBLAS.GxB_Matrix_apply_BinaryOp2nd(gbpointer(C), mask, accum, op, gbpointer(parent(A)), GBScalar(x), desc) return C end function apply!( - op, A::GBArray{T}, x; + op, A::GBArrayOrTranspose{T}, x; mask = nothing, accum = nothing, desc = nothing ) where {T} return apply!(op, A, A, x; mask, accum, desc) end function apply( - op, A::GBArray{T}, x; + op, A::GBArrayOrTranspose{T}, x; mask = nothing, accum = nothing, desc = nothing ) where {T} t = inferbinarytype(eltype(A), typeof(x), op) return apply!(op, similar(A, t), A, x; mask, accum, desc) end -function Base.map(f, A::GBArray{T}; mask = nothing, accum = nothing, desc = nothing) where {T} +function Base.map(f, A::GBArrayOrTranspose{T}; mask = nothing, accum = nothing, desc = nothing) where {T} apply(f, A; mask, accum, desc) end -function Base.map!(f, C::GBArray, A::GBArray{T}; mask = nothing, accum = nothing, desc = nothing) where {T} +function Base.map!(f, C::GBVecOrMat, A::GBArrayOrTranspose{T}; mask = nothing, accum = nothing, desc = nothing) where {T} apply!(f, C, A; mask, accum, desc) end -function Base.map!(f, A::GBArray{T}; mask = nothing, accum = nothing, desc = nothing) where {T} +function Base.map!(f, A::GBArrayOrTranspose{T}; mask = nothing, accum = nothing, desc = nothing) where {T} apply!(f, C, A; mask, accum, desc) end -Base.:*(x::V, u::GBArray{T}; mask = nothing, accum = nothing, desc = nothing) where {T, V<:Union{<:valid_union, T}} = +Base.:*(x::V, u::GBArrayOrTranspose{T}; mask = nothing, accum = nothing, desc = nothing) where {T, V<:Union{<:valid_union, T}} = apply(*, x, u; mask, accum, desc) -Base.:*(u::GBArray{T}, x::V; mask = nothing, accum = nothing, desc = nothing) where {T, V<:Union{<:valid_union, T}} = +Base.:*(u::GBArrayOrTranspose{T}, x::V; mask = nothing, accum = nothing, desc = nothing) where {T, V<:Union{<:valid_union, T}} = apply(*, u, x; mask, accum, desc) -Base.:-(u::GBArray) = apply(-, u) +Base.:-(u::GBArrayOrTranspose) = apply(-, u) """ - mask!(C::GBArray, A::GBArray, mask::GBArray) + mask!(C::GBArrayOrTranspose, A::GBArrayOrTranspose, mask::GBVecOrMat) Apply a mask to matrix `A`, storing the results in C. """ -function mask!(C::GBArray, A::GBArray, mask::GBArray; structural = false, complement = false) +function mask!(C::GBVecOrMat, A::GBArrayOrTranspose, mask::GBVecOrMat; structural = false, complement = false) desc = Descriptor() structural && (desc.structural_mask=true) complement && (desc.complement_mask=true) @@ -133,15 +133,15 @@ function mask!(C::GBArray, A::GBArray, mask::GBArray; structural = false, comple return C end -function mask!(A::GBArray, mask::GBArray; structural = false, complement = false) +function mask!(A::GBArrayOrTranspose, mask::GBVecOrMat; structural = false, complement = false) mask!(A, A, mask; structural, complement) end """ - mask(A::GBArray, mask::GBArray) + mask(A::GBArrayOrTranspose, mask::GBVecOrMat) Apply a mask to matrix `A`. """ -function mask(A::GBArray, mask::GBArray; structural = false, complement = false) +function mask(A::GBArrayOrTranspose, mask::GBVecOrMat; structural = false, complement = false) return mask!(similar(A), A, mask; structural, complement) end diff --git a/src/operations/mul.jl b/src/operations/mul.jl index c39b5d6d..78a78542 100644 --- a/src/operations/mul.jl +++ b/src/operations/mul.jl @@ -4,8 +4,8 @@ function LinearAlgebra.mul!( C::GBVecOrMat, - A::GBArray, - B::GBArray, + A::GBArrayOrTranspose, + B::GBArrayOrTranspose, op = (+, *); mask = nothing, accum = nothing, @@ -16,7 +16,7 @@ function LinearAlgebra.mul!( size(A, 2) == size(B, 1) || throw(DimensionMismatch("size(A, 2) != size(B, 1)")) size(A, 1) == size(C, 1) || throw(DimensionMismatch("size(A, 1) != size(C, 1)")) size(B, 2) == size(C, 2) || throw(DimensionMismatch("size(B, 2) != size(C, 2)")) - op = Semiring(op)(eltype(A), eltype(B)) + op = semiring(op, eltype(A), eltype(B)) accum = getaccum(accum, eltype(C)) op isa TypedSemiring || throw(ArgumentError("$op is not a valid TypedSemiring")) @wraperror LibGraphBLAS.GrB_mxm(gbpointer(C), mask, accum, op, gbpointer(parent(A)), gbpointer(parent(B)), desc) @@ -26,7 +26,7 @@ end function LinearAlgebra.mul!( C::GBVecOrMat, A::VecOrMat, - B::GBArray, + B::GBArrayOrTranspose, op = (+, *); mask = nothing, accum = nothing, @@ -40,7 +40,7 @@ end function LinearAlgebra.mul!( C::GBVecOrMat, - A::GBArray, + A::GBArrayOrTranspose, B::VecOrMat, op = (+, *); mask = nothing, @@ -54,7 +54,7 @@ function LinearAlgebra.mul!( end """ - *(A::GBArray, B::GBArray, op=(+,*); kwargs...)::GBArray + *(A::GBArrayOrTranspose, B::GBArrayOrTranspose, op=(+,*); kwargs...)::GBArrayOrTranspose Multiply two `GBArray`s `A` and `B` using a semiring, which defaults to the arithmetic semiring `+.*`. @@ -65,11 +65,11 @@ The mutating form, `mul!(C, A, B, op; kwargs...)` is identical except it stores The operator syntax `A * B` can be used when the default semiring is desired, and `*(max, +)(A, B)` can be used otherwise. # Arguments -- `A, B::GBArray`: A GBVector or GBMatrix, possibly transposed. +- `A, B::GBArrayOrTranspose`: A GBVector or GBMatrix, possibly transposed. - `op::Union{Tuple{Function, Function}, AbstractSemiring}`: the semiring used for matrix multiplication. May be passed as a tuple of functions, or an `AbstractSemiring` found in the `Semirings` submodule. # Keywords - `mask::Union{Nothing, GBArray} = nothing`: optional mask which determines the output pattern. -- `accum::Union{Nothing, Function, AbstractBinaryOp} = nothing`: optional binary accumulator +- `accum::Union{Nothing, Function} = nothing`: optional binary accumulator operation such that `C[i,j] = accum(C[i,j], T[i,j])` where T is the result of this function before accum is applied. - `desc::Union{Nothing, Descriptor}` @@ -78,8 +78,8 @@ The operator syntax `A * B` can be used when the default semiring is desired, an if a type specific semiring is provided. """ function Base.:*( - A::GBArray, - B::GBArray, + A::GBArrayOrTranspose, + B::GBArrayOrTranspose, op = (+, *); mask = nothing, accum = nothing, @@ -87,9 +87,9 @@ function Base.:*( ) t = inferbinarytype(eltype(A), eltype(B), op) fill = _promotefill(parent(A).fill, parent(B).fill) - if A isa GBMatOrTranspose && B isa AbstractGBVector + if A isa GBMatrixOrTranspose && B isa AbstractGBVector C = similar(A, t, size(A, 1); fill) - elseif A isa GBVector && B isa GBMatOrTranspose + elseif A isa GBVector && B isa GBMatrixOrTranspose C = similar(A, t, size(B, 2); fill) elseif A isa Transpose{<:Any, <:GBVector} && B isa GBVector C = similar(A, t, 1; fill) @@ -102,7 +102,7 @@ end function Base.:*( A::VecOrMat, - B::GBArray, + B::GBArrayOrTranspose, op = (+, *); mask = nothing, accum = nothing, @@ -115,7 +115,7 @@ function Base.:*( end function Base.:*( - A::GBArray, + A::GBArrayOrTranspose, B::VecOrMat, op = (+, *); mask = nothing, @@ -149,14 +149,20 @@ function Base.:*( return *(A, B, (+, *); mask, accum, desc) end -function Base.:*((⊕)::Function, (⊗)::Function) - return function(A::GBArray, B::GBArray; mask=nothing, accum=nothing, desc=nothing) +function Base.:*((⊕)::Union{<:Base.Callable, Monoid}, (⊗)::Function) + return function(A::GBArrayOrTranspose, B::GBArrayOrTranspose; mask=nothing, accum=nothing, desc=nothing) *(A, B, (⊕, ⊗); mask, accum, desc) end end -function Base.:*(rig::AbstractSemiring) - return function(A::GBArray, B::GBArray; mask=nothing, accum=nothing, desc=nothing) +function Base.:*((⊕)::Monoid, (⊗)::Function) + return function(A::GBArrayOrTranspose, B::GBArrayOrTranspose; mask=nothing, accum=nothing, desc=nothing) + *(A, B, (⊕, ⊗); mask, accum, desc) + end +end + +function Base.:*(rig::TypedSemiring) + return function(A::GBArrayOrTranspose, B::GBArrayOrTranspose; mask=nothing, accum=nothing, desc=nothing) *(A, B, rig; mask, accum, desc) end end diff --git a/src/operations/operationutils.jl b/src/operations/operationutils.jl index 9c1f72a9..bb6829e7 100644 --- a/src/operations/operationutils.jl +++ b/src/operations/operationutils.jl @@ -1,22 +1,20 @@ getaccum(::Nothing, t) = C_NULL getaccum(::Ptr{Nothing}, t) = C_NULL -getaccum(op::Function, t) = BinaryOp(op)(t, t) -getaccum(op::BinaryOp, t) = op(t, t) -getaccum(op::Function, tleft, tright) = BinaryOp(op)(tleft, tright) -getaccum(op::BinaryOp, tleft, tright) = op(tleft, tright) -getaccum(op::TypedBinaryOperator, tleft, tright=tleft) = op - -inferunarytype(::Type{T}, op::AbstractUnaryOp) where {T} = ztype(op(T)) -inferunarytype(::Type{T}, op) where {T} = inferunarytype(T, UnaryOp(op)) -inferunarytype(::Type{X}, op::TypedUnaryOperator{F, X, Z}) where {F, X, Z} = ztype(op) - -inferbinarytype(::Type{T}, ::Type{U}, op::AbstractBinaryOp) where {T, U} = ztype(op(T, U)) -inferbinarytype(::Type{T}, ::Type{U}, op) where {T, U} = inferbinarytype(T, U, BinaryOp(op)) -inferbinarytype(::Type{T}, ::Type{U}, op::AbstractMonoid) where {T, U} = inferbinarytype(T, U, op.binaryop) +getaccum(op::Function, t) = binaryop(op, t, t) +getaccum(op::Function, tleft, tright) = binaryop(op, tleft, tright) +getaccum(op::TypedBinaryOperator, x...) = op + +inferunarytype(::Type{T}, f::F) where {T, F<:Base.Callable} = Base._return_type(f, Tuple{T}) +inferunarytype(::Type{X}, op::TypedUnaryOperator) where X = ztype(op) + +inferbinarytype(::Type{T}, ::Type{U}, f::F) where {T, U, F<:Base.Callable} = Base._return_type(f, Tuple{T, U}) +# Overload for `first`, which will give Vector{T} normally: +inferbinarytype(::Type{T}, ::Type{U}, f::typeof(first)) where {T, U} = T +inferbinarytype(::Type{T}, ::Type{U}, op::AbstractMonoid) where {T, U} = inferbinarytype(T, U, op.binaryop.fn) #semirings are technically binary so we'll just overload that -inferbinarytype(::Type{T}, ::Type{U}, op::Tuple) where {T, U} = inferbinarytype(T, U, Semiring(op)) -inferbinarytype(::Type{T}, ::Type{U}, op::AbstractSemiring) where {T, U} = inferbinarytype(T, U, op.mulop) +inferbinarytype(::Type{T}, ::Type{U}, op::Tuple) where {T, U} = inferbinarytype(T, U, semiring(op, T, U)) +inferbinarytype(::Type{T}, ::Type{U}, op::TypedSemiring) where {T, U} = inferbinarytype(T, U, op.mulop) inferbinarytype(::Type{X}, ::Type{Y}, op::TypedBinaryOperator{F, X, Y, Z}) where {F, X, Y, Z} = ztype(op) inferbinarytype(::Type{X}, ::Type{X}, op::TypedMonoid{F, X, Z}) where {F, X, Z} = ztype(op) diff --git a/src/operations/reduce.jl b/src/operations/reduce.jl index 006f3d78..b6fd5695 100644 --- a/src/operations/reduce.jl +++ b/src/operations/reduce.jl @@ -1,10 +1,10 @@ function reduce!( - op, w::AbstractGBVector, A::GBMatOrTranspose; + op, w::AbstractGBVector, A::GBArrayOrTranspose; mask = nothing, accum = nothing, desc = nothing ) mask, accum = _handlenothings(mask, accum) desc = _handledescriptor(desc; in1=A) - op = Monoid(op)(eltype(w)) + op = typedmonoid(op, eltype(w)) accum = getaccum(accum, eltype(w)) @wraperror LibGraphBLAS.GrB_Matrix_reduce_Monoid( Ptr{LibGraphBLAS.GrB_Vector}(gbpointer(w)), mask, accum, op, gbpointer(parent(A)), desc @@ -14,7 +14,7 @@ end function Base.reduce( op, - A::GBArray; + A::GBArrayOrTranspose; dims = :, typeout = nothing, init = nothing, @@ -22,31 +22,32 @@ function Base.reduce( accum = nothing, desc = nothing ) + desc = _handledescriptor(desc; in1=A) mask, accum = _handlenothings(mask, accum) if typeout === nothing typeout = eltype(A) end - if dims == 2 && !(A isa GBVecOrTranspose) + if dims == 2 && !(A isa GBVectorOrTranspose) w = similar(A, typeout, size(A, 1)) reduce!(op, w, A; desc, accum, mask) return w - elseif dims == 1 && !(A isa GBVecOrTranspose) + elseif dims == 1 && !(A isa GBVectorOrTranspose) desc.transpose_input1 = true w = similar(A, typeout, size(A, 2)) reduce!(op, w, A; desc, accum, mask) return w - elseif dims == (1,2) || dims == Colon() || A isa GBVecOrTranspose + elseif dims == (1,2) || dims == Colon() || A isa GBVectorOrTranspose if init === nothing c = GBScalar{typeout}() - typec = typeout else - GBScalar{typeout}(init) - typec = typeof(init) + c = GBScalar{typeout}(init) + end + op = typedmonoid(op, typeout) + if nnz(c) == 1 && accum == C_NULL + accum = binaryop(op) end - op = Monoid(op)(typec) - desc = _handledescriptor(desc; in1=A) - accum = getaccum(accum, typec) + accum = getaccum(accum, typeout) @wraperror LibGraphBLAS.GrB_Matrix_reduce_Monoid_Scalar(c, accum, op, gbpointer(parent(A)), desc) return c[] end @@ -59,24 +60,24 @@ end Reduce `A` along dimensions of A with monoid `op`. # Arguments -- `op`: the reducer. This must map to an AbstractMonoid, not an AbstractBinaryOp. -- `A::GBArray`: `GBVector` or optionally transposed `GBMatrix`. +- `op`: the reducer. This must map to an AbstractMonoid, not a binary op. +- `A::GBArrayOrTranspose`: `GBVector` or optionally transposed `GBMatrix`. - `dims = :`: Optional dimensions for GBMatrix, may be `1`, `2`, or `:`. # Keywords - `typeout`: Optional output type specification. Defaults to `eltype(A)`. - `init`: Optional initial value. - `mask::Union{Nothing, GBMatrix} = nothing`: optional mask. -- `accum::Union{Nothing, AbstractBinaryOp} = nothing`: binary accumulator operation +- `accum::Union{Nothing} = nothing`: binary accumulator operation where `C[i,j] = accum(C[i,j], T[i,j])` where T is the result of this function before accum is applied. - `desc = nothing` """ reduce -Base.maximum(A::GBArray) = reduce(max, A) -Base.maximum(f::Function, A::GBArray) = reduce(max, map(f, A)) +Base.maximum(A::GBArrayOrTranspose) = reduce(max, A) +Base.maximum(f::Function, A::GBArrayOrTranspose) = reduce(max, map(f, A)) -Base.minimum(A::GBArray) = reduce(min, A) -Base.minimum(f::Function, A::GBArray) = reduce(min, map(f, A)) +Base.minimum(A::GBArrayOrTranspose) = reduce(min, A) +Base.minimum(f::Function, A::GBArrayOrTranspose) = reduce(min, map(f, A)) -Base.sum(A::GBArray; dims=:) = reduce(+, A; dims) +Base.sum(A::GBArrayOrTranspose; dims=:) = reduce(+, A; dims) diff --git a/src/operations/select.jl b/src/operations/select.jl index 5a4b0c51..861059b7 100644 --- a/src/operations/select.jl +++ b/src/operations/select.jl @@ -4,7 +4,7 @@ function select!( op, C::GBVecOrMat, - A::GBArray, + A::GBArrayOrTranspose, thunk = nothing; mask = nothing, accum = nothing, @@ -22,13 +22,13 @@ function select!( return C end -function select!(op, A::GBArray, thunk = nothing; mask = nothing, accum = nothing, desc = nothing) +function select!(op, A::GBArrayOrTranspose, thunk = nothing; mask = nothing, accum = nothing, desc = nothing) return select!(op, A, A, thunk; mask, accum, desc) end """ - select(op::Union{Function, SelectUnion}, A::GBArray; kwargs...)::GBArray - select(op::Union{Function, SelectUnion}, A::GBArray, thunk; kwargs...)::GBArray + select(op::Union{Function, SelectUnion}, A::GBArrayOrTranspose; kwargs...)::GBArrayOrTranspose + select(op::Union{Function, SelectUnion}, A::GBArrayOrTranspose, thunk; kwargs...)::GBArrayOrTranspose Return a `GBArray` whose elements satisfy the predicate defined by `op`. Some SelectOps or functions may require an additional argument `thunk`, for use in @@ -37,13 +37,13 @@ Some SelectOps or functions may require an additional argument `thunk`, for use # Arguments - `op::Union{Function, SelectUnion}`: A select operator from the SelectOps submodule. -- `A::GBArray` +- `A::GBArrayOrTranspose` - `thunk::Union{GBScalar, nothing, valid_union}`: Optional value used to evaluate `op`. # Keywords - `mask::Union{Nothing, GBMatrix} = nothing`: optional mask which determines the output pattern. -- `accum::Union{Nothing, AbstractBinaryOp} = nothing`: optional binary accumulator +- `accum::Union{Nothing} = nothing`: optional binary accumulator operation where `C[i,j] = accum(C[i,j], T[i,j])` where T is the result of this function before accum is applied. - `desc = nothing` @@ -52,7 +52,7 @@ Some SelectOps or functions may require an additional argument `thunk`, for use """ function select( op, - A::GBArray, + A::GBArrayOrTranspose, thunk = nothing; mask = nothing, accum = nothing, @@ -65,7 +65,7 @@ function select( return C end -LinearAlgebra.tril(A::GBArray, k::Integer = 0) = select(tril, A, k) -LinearAlgebra.triu(A::GBArray, k::Integer = 0) = select(triu, A, k) -SparseArrays.dropzeros(A::GBArray) = select(nonzeros, A) -SparseArrays.dropzeros!(A::GBArray) = select!(nonzeros, A) \ No newline at end of file +LinearAlgebra.tril(A::GBArrayOrTranspose, k::Integer = 0) = select(tril, A, k) +LinearAlgebra.triu(A::GBArrayOrTranspose, k::Integer = 0) = select(triu, A, k) +SparseArrays.dropzeros(A::GBArrayOrTranspose) = select(nonzeros, A) +SparseArrays.dropzeros!(A::GBArrayOrTranspose) = select!(nonzeros, A) \ No newline at end of file diff --git a/src/operations/sort.jl b/src/operations/sort.jl index febd0474..1602162d 100644 --- a/src/operations/sort.jl +++ b/src/operations/sort.jl @@ -1,20 +1,20 @@ # inplace sorts function Base.sort!( - C::Union{GBArray, Nothing}, #output matrix, possibly aliased with A, + C::Union{GBVecOrMat, Nothing}, #output matrix, possibly aliased with A, #or Nothing if just interested in P - P::Union{GBArray, Nothing}, # permutation matrix, possibly Nothing if just interested in C - A::GBArray; #input matrix, possibly aliased with C. + P::Union{GBVecOrMat, Nothing}, # permutation matrix, possibly Nothing if just interested in C + A::GBArrayOrTranspose; #input matrix, possibly aliased with C. dims = nothing, lt = <, desc = nothing # We only support a limited set of the keywords for now. # Missing by, rev, order, alg ) - A isa GBMatOrTranspose && dims === nothing && throw(ArgumentError("dims must be either 1 (sort columns) or 2 (sort rows) for matrix arguments.")) + A isa GBMatrixOrTranspose && dims === nothing && throw(ArgumentError("dims must be either 1 (sort columns) or 2 (sort rows) for matrix arguments.")) A isa GBVector && (dims = 1) C, P = _handlenothings(C, P) C == C_NULL && P == C_NULL && throw(ArgumentError("One (or both) of C and P must not be nothing.")) - op = BinaryOp(lt)(eltype(A)) + op = binaryop(lt, eltype(A)) if dims == 1 transpose = true elseif dims == 2 @@ -28,25 +28,25 @@ function Base.sort!( return C end -function Base.sort!(A::GBArray; dims = nothing, lt = <, desc = nothing) +function Base.sort!(A::GBVecOrMat; dims = nothing, lt = <, desc = nothing) return sort!(A, nothing, A; dims, lt, desc) end -function Base.sort!(C::GBArray, A::GBArray; dims, lt = <, desc = nothing) +function Base.sort!(C::GBVecOrMat, A::GBArrayOrTranspose; dims, lt = <, desc = nothing) return sort!(C, nothing, A; dims, lt, desc) end -function Base.sortperm!(P::GBArray, A::GBArray; dims = nothing, lt = <, desc = nothing) +function Base.sortperm!(P::GBVecOrMat, A::GBArrayOrTranspose; dims = nothing, lt = <, desc = nothing) sort!(nothing, P, A; dims, lt, desc) return P end -function Base.sort(A::GBArray; dims = nothing, lt = <, desc = nothing) +function Base.sort(A::GBArrayOrTranspose; dims = nothing, lt = <, desc = nothing) C = similar(A) return sort!(C, A; dims, lt, desc) end -function Base.sortperm(A::GBArray; dims = nothing, lt = <, desc = nothing) +function Base.sortperm(A::GBArrayOrTranspose; dims = nothing, lt = <, desc = nothing) P = similar(A, Int64) return sortperm!(P, A; dims, lt, desc) end diff --git a/src/operations/transpose.jl b/src/operations/transpose.jl index 54112bb7..ae979743 100644 --- a/src/operations/transpose.jl +++ b/src/operations/transpose.jl @@ -10,12 +10,12 @@ Eagerly evaluated matrix transpose, storing the output in `C`. # Keywords - `mask::Union{Nothing, GBMatrix} = nothing`: optional mask. -- `accum::Union{Nothing, AbstractBinaryOp} = nothing`: binary accumulator operation +- `accum::Union{Nothing} = nothing`: binary accumulator operation where `C[i,j] = accum(C[i,j], T[i,j])` where T is the result of this function before accum is applied. - `desc::Union{Nothing, Descriptor} = DEFAULTDESC` """ function gbtranspose!( - C::AbstractGBArray, A::GBArray; + C::AbstractGBArray, A::GBArrayOrTranspose; mask = nothing, accum = nothing, desc = nothing ) mask, accum = _handlenothings(mask, accum) @@ -31,20 +31,20 @@ Eagerly evaluated matrix transpose which returns the transposed matrix. # Keywords - `mask::Union{Nothing, GBMatrix} = nothing`: optional mask. -- `accum::Union{Nothing, AbstractBinaryOp} = nothing`: binary accumulator operation +- `accum::Union{Nothing} = nothing`: binary accumulator operation where `C[i,j] = accum(C[i,j], T[i,j])` where T is the result of this function before accum is applied. - `desc::Union{Nothing, Descriptor} = nothing` # Returns - `C::GBMatrix`: output matrix. """ -function gbtranspose(A::GBArray; mask = nothing, accum = nothing, desc = nothing) +function gbtranspose(A::GBArrayOrTranspose; mask = nothing, accum = nothing, desc = nothing) C = similar(A, size(A,2), size(A, 1)) gbtranspose!(C, A; mask, accum, desc) return C end -function LinearAlgebra.transpose(A::GBArray) +function LinearAlgebra.transpose(A::AbstractGBArray) return Transpose(A) end @@ -66,4 +66,5 @@ end LinearAlgebra.adjoint(A::GBVecOrMat) = transpose(A) #arrrrgh, type piracy. +# TODO: avoid this if possible LinearAlgebra.transpose(::Nothing) = nothing diff --git a/src/operators/binaryops.jl b/src/operators/binaryops.jl index 977934f1..d3eadb7c 100644 --- a/src/operators/binaryops.jl +++ b/src/operators/binaryops.jl @@ -1,45 +1,47 @@ module BinaryOps import ..SuiteSparseGraphBLAS -using ..SuiteSparseGraphBLAS: isGxB, isGrB, TypedBinaryOperator, AbstractBinaryOp, GBType, +using ..SuiteSparseGraphBLAS: isGxB, isGrB, TypedBinaryOperator, GBType, valid_vec, juliaop, gbtype, symtotype, Itypes, Ftypes, Ztypes, FZtypes, Rtypes, optype, Ntypes, Ttypes, suffix, valid_union using ..LibGraphBLAS -export BinaryOp, @binop +export BinaryOp, binaryop export second, rminus, iseq, isne, isgt, islt, isge, isle, ∨, ∧, lxor, xnor, fmod, bxnor, bget, bset, bclr, firsti0, firsti, firstj0, firstj, secondi0, secondi, secondj0, secondj, pair -struct BinaryOp{F} <: AbstractBinaryOp - juliaop::F +const BINARYOPS = IdDict{Tuple{<:Base.Callable, DataType, DataType}, TypedBinaryOperator}() + +function fallback_binaryop( + f::F, ::Type{X}, ::Type{Y} +) where {F<:Base.Callable, X, Y} + return get!(BINARYOPS, (f, X, Y)) do + TypedBinaryOperator(f, X, Y) + end end -SuiteSparseGraphBLAS.juliaop(op::BinaryOp) = op.juliaop -(op::BinaryOp)(T) = op(T, T) -BinaryOp(op::TypedBinaryOperator) = op +# If we have the same type we know we must fallback, +# more specific methods will be captured by dispatch. +binaryop(f::F, ::Type{X}, ::Type{X}) where {F<:Base.Callable, X} = fallback_binaryop(f, X, X) -# TODO, clean up this function, it allocates typedop and is otherwise perhaps a little slow. -function (op::BinaryOp)(::Type{T}, ::Type{U}; cont = true) where {T, U} #fallback - typedop = nothing - if T ∈ valid_vec && U ∈ valid_vec && cont # if both argtypes are in builtins *and* we're at 0th recursion: - promoted = optype(T, U) - if promoted ∈ valid_vec # check if we can promote within the builtins - typedop = try - op(promoted, promoted; cont=false) # try to find a builtin BinaryOp that matches. If we recurse back into this function we know there's no builtin, so cont=false - catch - nothing - end - end - end - typedop !== nothing && return typedop - resulttype = Base._return_type(op.juliaop, Tuple{T, U}) - if resulttype <: Tuple - throw(ArgumentError("Inferred a tuple return type for function $(string(op.juliaop)) on type $T.")) +function binaryop( + f::F, ::Type{X}, ::Type{Y} +) where {F<:Base.Callable, X, Y} + P = promote_type(X, Y) + if isconcretetype(P) + return binaryop(f, P, P) + else + return fallback_binaryop(f, X, Y) end - typedop = TypedBinaryOperator(op.juliaop, T, U, resulttype) - return typedop end +binaryop(f, type) = binaryop(f, type, type) +binaryop(op::TypedBinaryOperator, x...) = op + +SuiteSparseGraphBLAS.juliaop(op::TypedBinaryOperator) = op.fn + +# TODO, clean up this function, it allocates typedop and is otherwise perhaps a little slow. + function typedbinopconstexpr(jlfunc, builtin, namestr, xtype, ytype, outtype) # Complex ops must always be GxB prefixed if (xtype ∈ Ztypes || ytype ∈ Ztypes || outtype ∈ Ztypes) && isGrB(namestr) @@ -66,9 +68,9 @@ function typedbinopconstexpr(jlfunc, builtin, namestr, xtype, ytype, outtype) constquote = :(const $(esc(namesym)) = TypedBinaryOperator($(esc(jlfunc)), $(esc(xsym)), $(esc(ysym)), $(esc(outsym)))) end dispatchquote = if xtype === :Any && ytype === :Any - :((::$(esc(:BinaryOp)){$(esc(:typeof))($(esc(jlfunc)))})(::Type, ::Type; cont=false) = $(esc(namesym))) + :($(esc(:(SuiteSparseGraphBLAS.BinaryOps.binaryop)))(::$(esc(:typeof))($(esc(jlfunc))), ::Type, ::Type) = $(esc(namesym))) else - :((::$(esc(:BinaryOp)){$(esc(:typeof))($(esc(jlfunc)))})(::Type{$xsym}, ::Type{$ysym}; cont=false) = $(esc(namesym))) + :($(esc(:(SuiteSparseGraphBLAS.BinaryOps.binaryop)))(::$(esc(:typeof))($(esc(jlfunc))), ::Type{$xsym}, ::Type{$ysym}) = $(esc(namesym))) end return quote $(constquote) @@ -132,32 +134,54 @@ macro binop(expr...) end # All types @binop first GrB_FIRST T=>T -@binop new second GrB_SECOND T=>T +second(x, y) = y +@binop second GrB_SECOND T=>T @binop any GxB_ANY T=>T # this doesn't match the semantics of Julia's any, but that may be ok... -@binop new pair GrB_ONEB T=>T # I prefer pair, but to keep up with the spec I'll match... + +pair(x::T, y::T) where T = one(T) +@binop pair GrB_ONEB T=>T # I prefer pair, but to keep up with the spec I'll match... @binop (+) GrB_PLUS T=>T @binop (-) GrB_MINUS T=>T -@binop new rminus GxB_RMINUS T=>T + +rminus(x, y) = y - x +@binop rminus GxB_RMINUS T=>T + @binop (*) GrB_TIMES T=>T @binop (/) GrB_DIV T=>T @binop (\) GxB_RDIV T=>T @binop (^) GxB_POW T=>T + +iseq(x::T, y::T) where T = T(x == y) @binop new iseq GxB_ISEQ T=>T + +isne(x::T, y::T) where T = T(x != y) @binop new isne GxB_ISNE T=>T # Real types @binop min GrB_MIN R=>R @binop max GrB_MAX R=>R -@binop new isgt GxB_ISGT R=>R -@binop new islt GxB_ISLT R=>R -@binop new isge GxB_ISGE R=>R -@binop new isle GxB_ISLE R=>R -@binop new (∨) GxB_LOR R=>R -@binop new (∧) GxB_LAND R=>R -(::BinaryOp{typeof(|)})(::Type{Bool}, ::Type{Bool}) = LOR_BOOL -(::BinaryOp{typeof(&)})(::Type{Bool}, ::Type{Bool}) = LAND_BOOL -@binop new lxor GxB_LXOR R=>R +isgt(x::T, y::T) where T = T(x > y) +@binop isgt GxB_ISGT R=>R +islt(x::T, y::T) where T = T(x < y) +@binop islt GxB_ISLT R=>R +isge(x::T, y::T) where T = T(x >= y) +@binop isge GxB_ISGE R=>R +isle(x::T, y::T) where T = T(x <= y) +@binop isle GxB_ISLE R=>R +function ∨(x::T, y::T) where T + return (x != zero(T)) || (y != zero(T)) +end +@binop (∨) GxB_LOR R=>R +function ∧(x::T, y::T) where T + return (x != zero(T)) && (y != zero(T)) +end +@binop (∧) GxB_LAND R=>R +binaryop(::typeof(|), ::Type{Bool}, ::Type{Bool}) = LOR_BOOL +binaryop(::typeof(&), ::Type{Bool}, ::Type{Bool}) = LAND_BOOL + +lxor(x::T, y::T) where T = xor((x != zero(T)), (y != zero(T))) +@binop lxor GxB_LXOR R=>R # T/R => Bool @binop (==) GrB_EQ T=>Bool @@ -168,12 +192,13 @@ end @binop (<=) GrB_LE R=>Bool # Bool=>Bool, most of which are covered above. -@binop new xnor GrB_LXNOR Bool=>Bool +xnor(x::T, y::T) where T = !(lxor(x, y)) +@binop xnor GrB_LXNOR Bool=>Bool @binop atan GxB_ATAN2 F=>F @binop hypot GxB_HYPOT F=>F -@binop new fmod GxB_FMOD F=>F +@binop mod GxB_FMOD F=>F @binop rem GxB_REMAINDER F=>F @binop ldexp GxB_LDEXP F=>F @binop copysign GxB_COPYSIGN F=>F @@ -183,31 +208,39 @@ end @binop (|) GrB_BOR I=>I @binop (&) GrB_BAND I=>I @binop (⊻) GrB_BXOR I=>I -@binop new bxnor GrB_BXNOR I=>I -(::BinaryOp{typeof(⊻)})(::Type{Bool}, ::Type{Bool}) = LXOR_BOOL +bxnor(x::T, y::T) where T = ~⊻(x, y) +@binop bxnor GrB_BXNOR I=>I +binaryop(::typeof(⊻), ::Type{Bool}, ::Type{Bool}) = LXOR_BOOL +# leaving these without any equivalent Julia functions +# probably should only operate on Ints anyway. @binop new bget GxB_BGET I=>I @binop new bset GxB_BSET I=>I @binop new bclr GxB_BCLR I=>I @binop (>>) GxB_BSHIFT (I, Int8)=>I -# Positionals +# Positionals with dummy functions for output type inference purposes +firsti0(x, y) = 0::Int64 @binop new firsti0 GxB_FIRSTI Any=>N +firsti1(x, y) = 1::Int64 @binop new firsti GxB_FIRSTI1 Any=>N +firstj0(x, y) = 0::Int64 @binop new firstj0 GxB_FIRSTJ Any=>N +firstj1(x, y) = 1::Int64 @binop new firstj GxB_FIRSTJ1 Any=>N +secondi0(x, y) = 0::Int64 @binop new secondi0 GxB_SECONDI Any=>N +secondi1(x, y) = 1::Int64 @binop new secondi GxB_SECONDI1 Any=>N +secondj0(x, y) = 0::Int64 @binop new secondj0 GxB_SECONDJ Any=>N +secondj1(x, y) = 1::Int64 @binop new secondj GxB_SECONDJ1 Any=>N end - -const BinaryUnion = Union{AbstractBinaryOp, TypedBinaryOperator} - ztype(::TypedBinaryOperator{F, X, Y, Z}) where {F, X, Y, Z} = Z xtype(::TypedBinaryOperator{F, X, Y, Z}) where {F, X, Y, Z} = X ytype(::TypedBinaryOperator{F, X, Y, Z}) where {F, X, Y, Z} = Y diff --git a/src/operators/monoids.jl b/src/operators/monoids.jl index b027f42e..34326d6f 100644 --- a/src/operators/monoids.jl +++ b/src/operators/monoids.jl @@ -1,17 +1,42 @@ module Monoids import ..SuiteSparseGraphBLAS -using ..SuiteSparseGraphBLAS: isGxB, isGrB, TypedMonoid, AbstractMonoid, GBType, - valid_vec, juliaop, gbtype, symtotype, Itypes, Ftypes, Ztypes, FZtypes, Rtypes, Ntypes, Ttypes, suffix, BinaryOps.BinaryOp, _builtinMonoid, BinaryOps.∨, BinaryOps.∧, BinaryOps.lxor, BinaryOps.xnor, BinaryOps.bxnor +using ..SuiteSparseGraphBLAS: BinaryOps, isGxB, isGrB, TypedMonoid, AbstractMonoid, GBType, + valid_vec, juliaop, gbtype, symtotype, Itypes, Ftypes, Ztypes, FZtypes, Rtypes, + Ntypes, Ttypes, suffix, BinaryOps.binaryop, _builtinMonoid, BinaryOps.∨, + BinaryOps.∧, BinaryOps.lxor, BinaryOps.xnor, BinaryOps.bxnor, valid_union using ..LibGraphBLAS -export Monoid, @monoid +export Monoid, typedmonoid, defaultmonoid -struct Monoid{F} <: AbstractMonoid - binaryop::BinaryOp{F} +struct Monoid{F, I, T} <: AbstractMonoid + fn::F + identity::I + terminal::T end -SuiteSparseGraphBLAS.juliaop(op::Monoid) = juliaop(op.binaryop) -Monoid(f::Function) = Monoid(BinaryOp(f)) -Monoid(op::TypedMonoid) = op + +const MONOIDS = IdDict{Tuple{<:Monoid, DataType}, TypedMonoid}() + +SuiteSparseGraphBLAS.juliaop(op::Monoid) = op.fn + +function Monoid(fn, identity) + return Monoid(fn, identity, nothing) +end + +function typedmonoid(m::Monoid{F, I, Term}, ::Type{T}) where {F<:Base.Callable, I, Term, T} + return (get!(MONOIDS, (m, T)) do + TypedMonoid(binaryop(m.fn, T), m.identity, m.terminal) + end) +end + +typedmonoid(op::TypedMonoid, x...) = op + +# We default to no available monoid. +defaultmonoid(f::F, ::Type{T}) where {F<:Base.Callable, T} = nothing + +# Use defaultmonoid when available. User should verify that this results in the correct monoid. +typedmonoid(f::F, ::Type{T}) where {F<:Base.Callable, T} = typedmonoid(defaultmonoid(f, T), T) + +BinaryOps.binaryop(op::TypedMonoid) = op.binaryop # Can't really do ephemeral monoid fallbacks... Need the identity and possibly terminal. @@ -26,20 +51,19 @@ function typedmonoidconstexpr(jlfunc, builtin, namestr, type, identity, term) else namestr = namestr * "_$(suffix(type))" end - if builtin - namesym = Symbol(namestr[5:end]) - else - namesym = Symbol(namestr) - end typesym = Symbol(type) - if builtin - constquote = :(const $(esc(namesym)) = _builtinMonoid($namestr, BinaryOp($(esc(jlfunc)))($(esc(typesym)), $(esc(typesym))), $(esc(identity)), $(esc(term)))) - else - constquote = :(const $(esc(namesym)) = TypedMonoid(BinaryOp($(esc(jlfunc)))($(esc(typesym)), $(esc(typesym))), $(esc(identity)), $(esc(term)))) - end return quote - $(constquote) - (::$(esc(:Monoid)){$(esc(:typeof))($(esc(jlfunc)))})(::Type{$typesym}) = $(esc(namesym)) + $(esc(:MONOIDS))[ + ($(esc(:Monoid))($(esc(jlfunc)), $(esc(identity)), $(esc(term))), $(esc(typesym))) + ] = + _builtinMonoid( + $namestr, + binaryop($(esc(jlfunc)), + $(esc(typesym)), + $(esc(typesym))), + $(esc(identity)), + $(esc(term)) + ) end end @@ -104,28 +128,63 @@ end # We link to the BinaryOp rather than the Julia functions, # because users will mostly be exposed to the higher level interface. +const PLUSMONOID = Monoid(+, zero) +defaultmonoid(::typeof(+), ::Type{<:Number}) = PLUSMONOID @monoid (+) GrB_PLUS nB id=>zero -(::Monoid{typeof(+)})(::Type{Bool}) = LOR_MONOID_BOOL +# (::Monoid{typeof(+)})(::Type{Bool}) = LOR_MONOID_BOOL + +const TIMESMONOID = Monoid(*, one, zero) +defaultmonoid(::typeof(*), ::Type{<:Integer}) = TIMESMONOID @monoid (*) GrB_TIMES I id=>one term=>zero -(::Monoid{typeof(*)})(::Type{Bool}) = LAND_MONOID_BOOL +# (::Monoid{typeof(*)})(::Type{Bool}) = LAND_MONOID_BOOL +const FLOATTIMESMONOID = Monoid(*, one) # float * monoid doesn't have a terminal. +defaultmonoid(::typeof(*), ::Type{<:Union{Complex, AbstractFloat}}) = FLOATTIMESMONOID @monoid (*) GrB_TIMES FZ id=>one -@monoid any GxB_ANY T id=>one term=>one # This is technically incorrect. The identity and terminal are *ANY* value in the domain. +# This is technically incorrect. The identity and terminal are *ANY* value in the domain. +# TODO: Users MAY NOT extend the any monoid, and this should be banned somehow. +const ANYMONOID = Monoid(any, one, one) +defaultmonoid(::typeof(any), ::Type{<:valid_union}) = ANYMONOID +@monoid any GxB_ANY T id=>one term=>one + +const MINMONOID = Monoid(min, typemax, typemin) +defaultmonoid(::typeof(min), ::Type{<:Real}) = MINMONOID @monoid min GrB_MIN R id=>typemax term=>typemin + +const MAXMONOID = Monoid(max, typemin, typemax) +defaultmonoid(::typeof(max), ::Type{<:Real}) = MAXMONOID @monoid max GrB_MAX R id=>typemin term=>typemax +const ORMONOID = Monoid(∨, false, true) +defaultmonoid(::typeof(∨), ::Type{Bool}) = ORMONOID @monoid (∨) GrB_LOR Bool id=>false term=>true +const ANDMONOID = Monoid(∧, true, false) +defaultmonoid(::typeof(∧), ::Type{Bool}) = ANDMONOID @monoid (∧) GrB_LAND Bool id=>true term=>false + +const XORMONOID = Monoid(lxor, false) +defaultmonoid(::typeof(lxor), ::Type{Bool}) = XORMONOID @monoid (lxor) GrB_LXOR Bool id=>false + +const EQMONOID = Monoid(==, true) +defaultmonoid(::typeof(==), ::Type{Bool}) = EQMONOID @monoid (==) GrB_LXNOR Bool id=>true + +const BORMONOID = Monoid(|, zero, typemax) +defaultmonoid(::typeof(|), ::Type{<:Unsigned}) = BORMONOID @monoid (|) GrB_BOR I id=>zero term=>typemax +const BANDMONOID = Monoid(&, typemax, zero) +defaultmonoid(::typeof(&), ::Type{<:Unsigned}) = BANDMONOID @monoid (&) GrB_BAND I id=>typemax term=>zero +const BXORMONOID = Monoid(⊻, zero) +defaultmonoid(::typeof(⊻), ::Type{<:Unsigned}) = BXORMONOID @monoid (⊻) GrB_BXOR I id=>zero +const BXNORMONOID = Monoid(bxnor, typemax) +defaultmonoid(::typeof(bxnor), ::Type{<:Unsigned}) = BXORMONOID @monoid bxnor GrB_BXNOR I id=>typemax end -const MonoidUnion = Union{AbstractMonoid, TypedMonoid} ztype(::TypedMonoid{F, X, T}) where {F, X, T} = X xtype(::TypedMonoid{F, X, T}) where {F, X, T} = X ytype(::TypedMonoid{F, X, T}) where {F, X, T} = X \ No newline at end of file diff --git a/src/operators/operatorutils.jl b/src/operators/operatorutils.jl index ea339852..21a14f26 100644 --- a/src/operators/operatorutils.jl +++ b/src/operators/operatorutils.jl @@ -15,7 +15,7 @@ function optype(atype, btype) return promote_type(atype, btype) end end -optype(::GBArray{T}, ::GBArray{U}) where {T, U} = optype(T, U) +optype(::GBArrayOrTranspose{T}, ::GBArrayOrTranspose{U}) where {T, U} = optype(T, U) const Utypes = (:UInt8, :UInt16, :UInt32, :UInt64) const Itypes = (:Int8, :Int16, :Int32, :Int64, :UInt8, :UInt16, :UInt32, :UInt64) diff --git a/src/operators/selectops.jl b/src/operators/selectops.jl index f90c61a3..65284341 100644 --- a/src/operators/selectops.jl +++ b/src/operators/selectops.jl @@ -28,7 +28,7 @@ Select the entries **not** on the `k`th diagonal of A. """ function offdiag end #I don't know of a function which does this already. const OFFDIAG = SelectOp("GxB_OFFDIAG", LibGraphBLAS.GxB_SelectOp(), offdiag) -offdiag(A::GBArray, k=0) = select(offdiag, A, k) +offdiag(A::GBArrayOrTranspose, k=0) = select(offdiag, A, k) const NONZERO = SelectOp("GxB_NONZERO", LibGraphBLAS.GxB_SelectOp(), nonzeros) const EQ_ZERO = SelectOp("GxB_EQ_ZERO", LibGraphBLAS.GxB_SelectOp(), nothing) diff --git a/src/operators/semirings.jl b/src/operators/semirings.jl index cc58fa9e..dd5a9c42 100644 --- a/src/operators/semirings.jl +++ b/src/operators/semirings.jl @@ -1,28 +1,35 @@ module Semirings import ..SuiteSparseGraphBLAS -using ..SuiteSparseGraphBLAS: isGxB, isGrB, TypedSemiring, AbstractSemiring, GBType, +using ..SuiteSparseGraphBLAS: isGxB, isGrB, TypedSemiring, GBType, valid_vec, juliaop, gbtype, symtotype, Itypes, Ftypes, Ztypes, FZtypes, - Rtypes, Ntypes, Ttypes, suffix, BinaryOps.BinaryOp, Monoids.Monoid, BinaryOps.second, BinaryOps.rminus, BinaryOps.pair, + Rtypes, Ntypes, Ttypes, suffix, BinaryOps.binaryop, Monoids.Monoid, BinaryOps.second, BinaryOps.rminus, BinaryOps.pair, BinaryOps.iseq, BinaryOps.isne, BinaryOps.isgt, BinaryOps.islt, BinaryOps.isge, BinaryOps.isle, BinaryOps.∨, - BinaryOps.∧, BinaryOps.lxor, BinaryOps.xnor, BinaryOps.fmod, BinaryOps.bxnor, BinaryOps.bget, BinaryOps.bset, + BinaryOps.∧, BinaryOps.lxor, BinaryOps.xnor, mod, BinaryOps.bxnor, BinaryOps.bget, BinaryOps.bset, BinaryOps.bclr, BinaryOps.firsti0, BinaryOps.firsti, BinaryOps.firstj0, BinaryOps.firstj, BinaryOps.secondi0, - BinaryOps.secondi, BinaryOps.secondj0, BinaryOps.secondj, xtype, ytype, ztype + BinaryOps.secondi, BinaryOps.secondj0, BinaryOps.secondj, xtype, ytype, ztype, + Monoids.typedmonoid, Monoids.defaultmonoid, valid_union using ..LibGraphBLAS -export Semiring, @rig +export semiring -struct Semiring{FM, FA} <: AbstractSemiring - addop::Monoid{FM} - mulop::BinaryOp{FA} -end -Semiring(addop::Function, mulop::Function) = Semiring(Monoid(addop),BinaryOp(mulop)) -Semiring(tup::Tuple{Function, Function}) = Semiring(tup...) -Semiring(op::TypedSemiring) = op -(rig::Semiring)(type) = rig(type, type) -function (rig::Semiring)(T, U) #fallback - mulop = rig.mulop(T, U) - Z = ztype(mulop) - TypedSemiring(rig.addop(Z), mulop) +const SEMIRINGS = IdDict{Tuple{<:Monoid, <:Base.Callable, DataType, DataType}, TypedSemiring}() + +function semiring( + m::Monoid{F, I, Term}, f::F2, ::Type{T}, ::Type{T2} +) where {F<:Base.Callable, F2<:Base.Callable, I, Term, T, T2} + return (get!(SEMIRINGS, (m, f, T, T2)) do + TypedSemiring(typedmonoid(m, T), binaryop(f, T, T2)) + end) end + +semiring(addop::Function, mulop::Function, ::Type{T}, ::Type{U}) where {T,U} = + semiring(defaultmonoid(addop, Base._return_type(mulop, Tuple{T, U})), mulop, T, U) +semiring(tup::Tuple{<:Union{Function, Monoid}, Function}, ::Type{T}, ::Type{U}) where {T,U} = semiring(tup..., T, U) +semiring(op::TypedSemiring, ::Type{T}, ::Type{U}) where {T,U} = op + +semiring(addop::Function, mulop::Function, ::Type{T}) where T = semiring(defaultmonoid(addop),mulop, T, T) +semiring(tup::Tuple{Function, Function}, ::Type{T}) where T = semiring(tup..., T, T) +semiring(op::TypedSemiring, ::Type{T}) where T = op + function typedrigconstexpr(addfunc, mulfunc, builtin, namestr, xtype, ytype, outtype) # Complex ops must always be GxB prefixed if (xtype ∈ Ztypes || ytype ∈ Ztypes || outtype ∈ Ztypes) && isGrB(namestr) @@ -45,12 +52,32 @@ function typedrigconstexpr(addfunc, mulfunc, builtin, namestr, xtype, ytype, out ysym = Symbol(ytype) outsym = Symbol(outtype) dispatchquote = if xtype === :Any && ytype === :Any - :((::$(esc(:Semiring)){$(esc(:typeof))($(esc(addfunc))), $(esc(:typeof))($(esc(mulfunc)))})(::Type, ::Type) = $(esc(namesym))) + quote + $(esc(:semiring))( + ::$(esc(:typeof))($(esc(addfunc))), + ::$(esc(:typeof))($(esc(mulfunc))), + ::Type{T}, + ::Type{T} + ) where {T<:valid_union}= $(esc(namesym)) + end else - :((::$(esc(:Semiring)){$(esc(:typeof))($(esc(addfunc))), $(esc(:typeof))($(esc(mulfunc)))})(::Type{$xsym}, ::Type{$ysym}) = $(esc(namesym))) + quote + $(esc(:semiring))( + ::$(esc(:typeof))($(esc(addfunc))), + ::$(esc(:typeof))($(esc(mulfunc))), + ::Type{$xsym}, ::Type{$ysym} + ) = $(esc(namesym)) + end end return quote - const $(esc(namesym)) = TypedSemiring($builtin, false, $namestr, LibGraphBLAS.GrB_Semiring(), Monoid($(esc(addfunc)))($(esc(outsym))), BinaryOp($(esc(mulfunc)))($(esc(xsym)), $(esc(ysym)))) + const $(esc(namesym)) = TypedSemiring( + $builtin, + false, + $namestr, + LibGraphBLAS.GrB_Semiring(), + typedmonoid(defaultmonoid($(esc(addfunc)), $(esc(outsym))), $(esc(outsym))), + binaryop($(esc(mulfunc)), $(esc(xsym)), $(esc(ysym))) + ) $(dispatchquote) end end diff --git a/src/operators/unaryops.jl b/src/operators/unaryops.jl index f1bce576..be786f36 100644 --- a/src/operators/unaryops.jl +++ b/src/operators/unaryops.jl @@ -1,32 +1,27 @@ module UnaryOps import ..SuiteSparseGraphBLAS -using ..SuiteSparseGraphBLAS: isGxB, isGrB, TypedUnaryOperator, AbstractUnaryOp, GBType, +using ..SuiteSparseGraphBLAS: isGxB, isGrB, TypedUnaryOperator, GBType, valid_vec, juliaop, gbtype, symtotype, Itypes, Ftypes, Ztypes, FZtypes, Rtypes, Ntypes, Ttypes, suffix using ..LibGraphBLAS -export UnaryOp, @unop +export unaryop, @unop export positioni, positionj, frexpx, frexpe using SpecialFunctions -struct UnaryOp{F} <: AbstractUnaryOp - juliaop::F +const UNARYOPS = IdDict{Tuple{<:Base.Callable, DataType}, TypedUnaryOperator}() + +function unaryop(f::F, ::Type{T}) where {F<:Base.Callable, T} + return get!(UNARYOPS, (f, T)) do + TypedUnaryOperator(f, T) + end end -UnaryOp(op::TypedUnaryOperator) = op +unaryop(op::TypedUnaryOperator, x...) = op -SuiteSparseGraphBLAS.juliaop(op::UnaryOp) = op.juliaop +SuiteSparseGraphBLAS.juliaop(op::TypedUnaryOperator) = op.fn -# fallback -# This is a fallback since creating these every time is incredibly costly. -function (op::UnaryOp)(::Type{T}) where {T} - resulttype = Base._return_type(op.juliaop, Tuple{T}) - if resulttype <: Tuple - throw(ArgumentError("Inferred a tuple return type for function $(string(op.juliaop)) on type $T.")) - end - return TypedUnaryOperator(op.juliaop, T, resulttype) -end function typedunopconstexpr(jlfunc, builtin, namestr, intype, outtype) # Complex ops must always be GxB prefixed if (intype ∈ Ztypes || outtype ∈ Ztypes) && isGrB(namestr) @@ -51,7 +46,7 @@ function typedunopconstexpr(jlfunc, builtin, namestr, intype, outtype) end return quote $(constquote) - (::$(esc(:UnaryOp)){$(esc(:typeof))($(esc(jlfunc)))})(::Type{$(esc(insym))}) = $(esc(namesym)) + $(esc(:(SuiteSparseGraphBLAS.UnaryOps.unaryop)))(::$(esc(:typeof))($(esc(jlfunc))), ::Type{$(esc(insym))}) = $(esc(namesym)) end end @@ -104,7 +99,7 @@ macro unop(expr...) end # all types -@unop identity GrB_IDENTITY Bool=>Bool +@unop identity GrB_IDENTITY T=>T @unop (-) GrB_AINV T=>T @unop inv GrB_MINV T=>T @unop one GxB_ONE T=>T @@ -115,8 +110,13 @@ end @unop (~) GrB_BNOT I=>I # positionals -@unop new positioni GxB_POSITIONI1 Any=>N -@unop new positionj GxB_POSITIONJ1 Any=>N +# dummy functions mostly for Base._return_type purposes. +# 1 is the most natural value regardless. +positioni(_) = 1::Int64 +positionj(_) = 1::Int64 +@unop positioni GxB_POSITIONI1 Any=>N +@unop positionj GxB_POSITIONJ1 Any=>N + #floats and complexes @unop sqrt GxB_SQRT FZ=>FZ @@ -153,12 +153,43 @@ end @unop erf GxB_ERF FZ=>FZ @unop erfc GxB_ERFC FZ=>FZ # julia has frexp which returns (x, exp). This is split in SS:GrB to frexpx = frexp[1]; frexpe = frexp[2]; -@unop new frexpx GxB_FREXPX FZ=>FZ -@unop new frexpe GxB_FREXPE FZ=>FZ +frexpx(x) = frexp[1] +frexpe(x) = frexp[2] +@unop frexpx GxB_FREXPX FZ=>FZ +@unop frexpe GxB_FREXPE FZ=>FZ @unop isinf GxB_ISINF FZ=>Bool @unop isnan GxB_ISNAN FZ=>Bool @unop isfinite GxB_ISFINITE FZ=>Bool +# manually create promotion overloads. +# Otherwise we will fallback to Julia functions which can be harmful. +for f ∈ [ + sqrt, log, exp, log10, log2, exp2, expm1, log1p, sin, cos, tan, + asin, acos, atan, sinh, cosh, tanh, asinh, acosh, atanh] + @eval begin + unaryop(::typeof($f), ::Type{UInt64}) = + unaryop($f, Float64) + unaryop(::typeof($f), ::Type{<:Union{UInt32, UInt16, UInt8}}) = + unaryop($f, Float32) + end +end + +# I think this list is correct. +# It should be those functions which can be safely promoted to float +# from Ints (including negatives). +# This might be overzealous, and should be just combined with the list above. +# I'd rather error on domain than create a bunch of NaNs. +for f ∈ [ + exp, exp2, expm1, sin, cos, tan, + atan, sinh, cosh, tanh, asinh] + @eval begin + unaryop(::typeof($f), ::Type{Int64}) = + unaryop($f, Float64) + unaryop(::typeof($f), ::Type{<:Union{Int32, Int16, Int8}}) = + unaryop($f, Float32) + end +end + # Complex functions @unop conj GxB_CONJ Z=>Z @unop real GxB_CREAL Z=>F @@ -166,8 +197,5 @@ end @unop angle GxB_CARG Z=>F end -const UnaryUnion = Union{AbstractUnaryOp, TypedUnaryOperator} - - ztype(::TypedUnaryOperator{F, I, O}) where {F, I, O} = O xtype(::TypedUnaryOperator{F, I, O}) where {F, I, O} = I diff --git a/src/scalar.jl b/src/scalar.jl index 6cf68539..a8e16cfb 100644 --- a/src/scalar.jl +++ b/src/scalar.jl @@ -12,6 +12,18 @@ function GBScalar(v::T) where {T} return x end +function GBScalar{T}(v::T) where {T} + x = GBScalar{T}() + x[] = v + return x +end + +function GBScalar{T}(v::U) where {T, U} + x = GBScalar{T}() + x[] = convert(T, v) + return x +end + # Some Base and basic SparseArrays/LinearAlgebra functions: ########################################################### Base.unsafe_convert(::Type{LibGraphBLAS.GxB_Scalar}, s::GBScalar) = s.p @@ -65,5 +77,5 @@ end function SparseArrays.nnz(v::GBScalar) n = Ref{LibGraphBLAS.GrB_Index}() @wraperror LibGraphBLAS.GrB_Scalar_nvals(n, v) - return n[] + return Int64(n[]) end diff --git a/src/types.jl b/src/types.jl index be77c215..337d95d4 100644 --- a/src/types.jl +++ b/src/types.jl @@ -26,6 +26,10 @@ function TypedUnaryOperator(fn::F, ::Type{X}, ::Type{Z}) where {F, X, Z} return TypedUnaryOperator{F, X, Z}(false, false, string(fn), LibGraphBLAS.GrB_UnaryOp(), fn) end +function TypedUnaryOperator(fn::F, ::Type{X}) where {F, X} + return TypedUnaryOperator(fn, X, Base._return_type(fn, Tuple{X})) +end + function Base.unsafe_convert(::Type{LibGraphBLAS.GrB_UnaryOp}, op::TypedUnaryOperator{F, X, Z}) where {F, X, Z} # We can lazily load the built-ins since they are already constants. # Could potentially do this with UDFs, but probably not worth the effort. @@ -71,6 +75,10 @@ function TypedBinaryOperator(fn::F, ::Type{X}, ::Type{Y}, ::Type{Z}) where {F, X return TypedBinaryOperator{F, X, Y, Z}(false, false, string(fn), LibGraphBLAS.GrB_BinaryOp(), fn) end +function TypedBinaryOperator(fn::F, ::Type{X}, ::Type{Y}) where {F, X, Y} + return TypedBinaryOperator(fn, X, Y, Base._return_type(fn, Tuple{X, Y})) +end + function (op::TypedBinaryOperator{F, X, Y, Z})(::Type{T1}, ::Type{T2}) where {F, X, Y, Z, T1, T2} return op end @@ -168,6 +176,7 @@ function Base.unsafe_convert(::Type{LibGraphBLAS.GrB_Monoid}, op::TypedMonoid{F, end end +# TODO: Determine TERMINAL::T, T ∈ Z always, or if T can be outside Z. function _builtinMonoid(typestr, binaryop::TypedBinaryOperator{F, Z, Z, Z}, identity, terminal::T) where {F, Z, T} return TypedMonoid(true, false, typestr, LibGraphBLAS.GrB_Monoid(), binaryop, identity, terminal) end @@ -176,6 +185,10 @@ function TypedMonoid(binop::TypedBinaryOperator{F, Z, Z, Z}, identity::Z, termin return TypedMonoid(false, false, string(binop.fn), LibGraphBLAS.GrB_Monoid(), binop, identity, terminal) end +function TypedMonoid(binop::TypedBinaryOperator{F, Z, Z, Z}, identity, terminal::T) where {F, Z, T} + return TypedMonoid(false, false, string(binop.fn), LibGraphBLAS.GrB_Monoid(), binop, convert(Z, identity), terminal) +end + #Enable use of functions for determining identity and terminal values. Could likely be pared down to 2 functions somehow. TypedMonoid(builtin, loaded, typestr, p, binaryop::TypedBinaryOperator{F, Z, Z, Z}, identity::Function, terminal::Function) where {F, Z} = TypedMonoid(builtin, loaded, typestr, p, binaryop, identity(Z), terminal(Z)) diff --git a/src/unpack.jl b/src/unpack.jl index be313f64..4bf33bb6 100644 --- a/src/unpack.jl +++ b/src/unpack.jl @@ -123,10 +123,10 @@ function _unpackcscmatrix!(A::AbstractGBArray{T}; desc = nothing, incrementindic colvector = finalizer(unsafe_wrap(Array, Ptr{Int64}(colptr[]), size(A, 2) + 1)) do x _jlfree(x) end - rowvector = finalizer(unsafe_wrap(Array, Ptr{Int64}(rowidx[]), rowidxsize[])) do x + rowvector = finalizer(unsafe_wrap(Array, Ptr{Int64}(rowidx[]), nnonzeros)) do x _jlfree(x) end - v = finalizer(unsafe_wrap(Array, Ptr{T}(values[]), nnz)) do x + v = finalizer(unsafe_wrap(Array, Ptr{T}(values[]), nnonzeros)) do x _jlfree(x) end x = (colvector, rowvector, v) diff --git a/test/gbarray.jl b/test/gbarray.jl index 3b712f43..8f123ae7 100644 --- a/test/gbarray.jl +++ b/test/gbarray.jl @@ -17,8 +17,8 @@ x = sprand(Int64, 100, 100, 0.05) m = GBMatrix(x) @test m[1, 2] === nothing - @test m[:, 2] == GBMatrix(x[:, 2]) - @test m[2, :] == copy(GBMatrix(x[2, :])') + @test m[:, 2] == GBVector(x[:, 2]) + @test m[2, :] == copy(GBVector(x[2, :])') @test m[:, :] == m @test m[1:2:5, 1:2] == GBMatrix(x[1:2:5, 1:2]) @test m[1:2:5, :] == GBMatrix(x[1:2:5, :]) diff --git a/test/ops.jl b/test/ops.jl index 04383d1f..f4c46fd2 100644 --- a/test/ops.jl +++ b/test/ops.jl @@ -1,16 +1,14 @@ @testset "Ops" begin @testset "UnaryOps" begin - @test SuiteSparseGraphBLAS.SuiteSparseGraphBLAS.juliaop(UnaryOp(sin)) === sin - @test UnaryOp(sin)(Float64).builtin + @test SuiteSparseGraphBLAS.SuiteSparseGraphBLAS.juliaop(unaryop(sin, Float64)) === sin + @test unaryop(sin, Float64).builtin #test the manual construction of TypedUnaryOperators X = GBVector(Float64[1,3,5]) f = (x) -> 1.3 # a random unary function. - op = UnaryOp(f) - @test SuiteSparseGraphBLAS.SuiteSparseGraphBLAS.juliaop(op) == f - typedop = op(Float64) + typedop = unaryop(f, Float64) @test !typedop.builtin @test !typedop.loaded - @test UnaryOp(typedop) == typedop + @test unaryop(typedop, Float64) == typedop @test map(typedop, X)[1] == 1.3 @test typedop.loaded @@ -24,22 +22,18 @@ end @testset "BinaryOps" begin # kinda vacuous tests here... - @test xtype(BinaryOp(+)(Float64)) == Float64 - @test ytype(BinaryOp(+)(Float64)) == Float64 - @test ztype(BinaryOp(+)(Float64)) == Float64 + @test xtype(SuiteSparseGraphBLAS.binaryop(+, Float64)) == Float64 + @test ytype(SuiteSparseGraphBLAS.binaryop(+, Float64)) == Float64 + @test ztype(SuiteSparseGraphBLAS.binaryop(+, Float64)) == Float64 end @testset "Monoids" begin - @test xtype(Monoid(+)(Float64)) == Float64 - @test ytype(Monoid(+)(Float64)) == Float64 - @test ztype(Monoid(+)(Float64)) == Float64 - - op = Monoid(+) - @test SuiteSparseGraphBLAS.juliaop(op) === + - @test Monoid(op(ComplexF64)) == op(ComplexF64) + @test xtype(SuiteSparseGraphBLAS.typedmonoid(+, Float64)) == Float64 + @test ytype(SuiteSparseGraphBLAS.typedmonoid(+, Float64)) == Float64 + @test ztype(SuiteSparseGraphBLAS.typedmonoid(+, Float64)) == Float64 end @testset "Semirings" begin - @test xtype(Semiring(+, *)(ComplexF64)) == ComplexF64 - @test ytype(Semiring(+, *)(ComplexF64)) == ComplexF64 - @test ztype(Semiring(+, *)(ComplexF64)) == ComplexF64 + @test xtype(semiring((+, *), ComplexF64)) == ComplexF64 + @test ytype(semiring((+, *), ComplexF64)) == ComplexF64 + @test ztype(semiring((+, *), ComplexF64)) == ComplexF64 end end \ No newline at end of file