diff --git a/src/SuiteSparseGraphBLAS.jl b/src/SuiteSparseGraphBLAS.jl index 5663bee8..2de5ccb2 100644 --- a/src/SuiteSparseGraphBLAS.jl +++ b/src/SuiteSparseGraphBLAS.jl @@ -55,6 +55,7 @@ include("operations/extract.jl") include("scalar.jl") include("vector.jl") include("matrix.jl") +include("abstractgbarray.jl") include("random.jl") # include("operations/operationutils.jl") @@ -71,7 +72,6 @@ include("operations/sort.jl") # include("print.jl") include("import.jl") -include("export.jl") include("pack.jl") include("unpack.jl") include("options.jl") @@ -88,7 +88,7 @@ include("chainrules/constructorrules.jl") #EXPERIMENTAL include("misc.jl") include("asjulia.jl") -include("spmgb/sparsemat.jl") +# include("spmgb/sparsemat.jl") include("mmread.jl") export SparseArrayCompat diff --git a/src/abstractgbarray.jl b/src/abstractgbarray.jl new file mode 100644 index 00000000..8fa9c461 --- /dev/null +++ b/src/abstractgbarray.jl @@ -0,0 +1,527 @@ +# AbstractGBArray functions: +function SparseArrays.nnz(A::AbsGBArrayOrTranspose) + 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 + +""" + clear!(v::GBVector) + clear!(A::GBMatrix) + +Clear all the entries from the GBArray. +Does not modify the type or dimensions. +""" +clear!(A::AbsGBArrayOrTranspose) = @wraperror LibGraphBLAS.GrB_Matrix_clear(gbpointer(parent(A))); return nothing + +# AbstractGBMatrix functions: +############################# + +function build(A::AbstractGBMatrix{T}, I::AbstractVector, J::AbstractVector, x::T) where {T} + nnz(A) == 0 || throw(OutputNotEmptyError("Cannot build matrix with existing elements")) + length(I) == length(J) || DimensionMismatch("I, J and X must have the same length") + x = GBScalar(x) + + @wraperror LibGraphBLAS.GxB_Matrix_build_Scalar( + gbpointer(A), + Vector{LibGraphBLAS.GrB_Index}(decrement!(I)), + Vector{LibGraphBLAS.GrB_Index}(decrement!(J)), + x, + length(I) + ) + increment!(I) + increment!(J) + return A +end + +function Base.size(A::AbstractGBMatrix) + nrows = Ref{LibGraphBLAS.GrB_Index}() + ncols = Ref{LibGraphBLAS.GrB_Index}() + @wraperror LibGraphBLAS.GrB_Matrix_nrows(nrows, gbpointer(A)) + @wraperror LibGraphBLAS.GrB_Matrix_ncols(ncols, gbpointer(A)) + return (Int64(nrows[]), Int64(ncols[])) +end + +function Base.deleteat!(A::AbstractGBMatrix, i, j) + @wraperror LibGraphBLAS.GrB_Matrix_removeElement(A, decrement!(i), decrement!(j)) + return A +end + +function Base.resize!(A::AbstractGBMatrix, nrows_new, ncols_new) + @wraperror LibGraphBLAS.GrB_Matrix_resize(gbpointer(A), nrows_new, ncols_new) + return A +end + +# Type dependent functions build, setindex, getindex, and findnz: +for T ∈ valid_vec + if T ∈ gxb_vec + prefix = :GxB + else + prefix = :GrB + end + # Build functions + func = Symbol(prefix, :_Matrix_build_, suffix(T)) + @eval begin + function build(A::AbstractGBMatrix{$T}, I::AbstractVector{<:Integer}, J::AbstractVector{<:Integer}, X::AbstractVector{$T}; + combine = + + ) + combine = BinaryOp(combine)($T) + I isa Vector || (I = collect(I)) + J isa Vector || (J = collect(J)) + X isa Vector || (X = collect(X)) + nnz(A) == 0 || throw(OutputNotEmptyError("Cannot build matrix with existing elements")) + length(X) == length(I) == length(J) || + DimensionMismatch("I, J and X must have the same length") + decrement!(I) + decrement!(J) + @wraperror LibGraphBLAS.$func( + gbpointer(A), + I, + J, + X, + length(X), + combine + ) + increment!(I) + increment!(J) + end + end + # Setindex functions + func = Symbol(prefix, :_Matrix_setElement_, suffix(T)) + @eval begin + function Base.setindex!(A::AbstractGBMatrix{$T}, x, i::Integer, j::Integer) + x = convert($T, x) + @wraperror LibGraphBLAS.$func(gbpointer(A), x, LibGraphBLAS.GrB_Index(decrement!(i)), LibGraphBLAS.GrB_Index(decrement!(j))) + return x + end + end + # Getindex functions + func = Symbol(prefix, :_Matrix_extractElement_, suffix(T)) + @eval begin + function Base.getindex(A::AbstractGBMatrix{$T}, i::Int, j::Int) + x = Ref{$T}() + result = LibGraphBLAS.$func(x, gbpointer(A), decrement!(i), decrement!(j)) + if result == LibGraphBLAS.GrB_SUCCESS + return x[] + elseif result == LibGraphBLAS.GrB_NO_VALUE + return A.fill + else + @wraperror result + end + end + # Fix ambiguity + function Base.getindex(A::Transpose{$T, <:AbstractGBMatrix{$T}}, i::Int, j::Int) + return getindex(parent(A), j, i) + end + end + # findnz functions + func = Symbol(prefix, :_Matrix_extractTuples_, suffix(T)) + @eval begin + function SparseArrays.findnz(A::AbstractGBMatrix{$T}) + nvals = Ref{LibGraphBLAS.GrB_Index}(nnz(A)) + I = Vector{LibGraphBLAS.GrB_Index}(undef, nvals[]) + J = Vector{LibGraphBLAS.GrB_Index}(undef, nvals[]) + X = Vector{$T}(undef, nvals[]) + wait(A) + @wraperror LibGraphBLAS.$func(I, J, X, nvals, gbpointer(A)) + nvals[] == length(I) == length(J) == length(X) || throw(DimensionMismatch("length(I) != length(X)")) + return increment!(I), increment!(J), X + end + function SparseArrays.nonzeros(A::AbstractGBMatrix{$T}) + nvals = Ref{LibGraphBLAS.GrB_Index}(nnz(A)) + X = Vector{$T}(undef, nvals[]) + wait(A) + @wraperror LibGraphBLAS.$func(C_NULL, C_NULL, X, nvals, gbpointer(A)) + nvals[] == length(X) || throw(DimensionMismatch("")) + return X + end + function SparseArrays.nonzeroinds(A::AbstractGBMatrix{$T}) + nvals = Ref{LibGraphBLAS.GrB_Index}(nnz(A)) + I = Vector{LibGraphBLAS.GrB_Index}(undef, nvals[]) + J = Vector{LibGraphBLAS.GrB_Index}(undef, nvals[]) + wait(A) + @wraperror LibGraphBLAS.$func(I, J, C_NULL, nvals, gbpointer(A)) + nvals[] == length(I) == length(J) || throw(DimensionMismatch("")) + return increment!(I), increment!(J) + end + end +end + +for T ∈ valid_vec + func = Symbol(:GxB_Matrix_subassign_, suffix(T)) + @eval begin + function _subassign(C::AbstractGBMatrix{$T}, x, I, ni, J, nj, mask, accum, desc) + @wraperror LibGraphBLAS.$func(gbpointer(C), mask, accum, x, I, ni, J, nj, desc) + return x + end + end + if T ∈ gxb_vec + prefix = :GxB + else + prefix = :GrB + end + func = Symbol(prefix, :_Matrix_assign_, suffix(T)) + @eval begin + function _assign(C::AbstractGBMatrix{$T}, x, I, ni, J, nj, mask, accum, desc) + @wraperror LibGraphBLAS.$func(C, mask, accum, x, I, ni, J, nj, desc) + return x + end + end + # TODO: Update when upstream. + # this is less than ideal. But required for isstored. + # a new version of graphBLAS will replace this with Matrix_extractElement_Structural + func = Symbol(prefix, :_Matrix_extractElement_, suffix(T)) + @eval begin + function Base.isstored(A::AbstractGBMatrix{$T}, i::Int, j::Int) + result = LibGraphBLAS.$func(Ref{$T}(), gbpointer(A), decrement!(i), decrement!(j)) + if result == LibGraphBLAS.GrB_SUCCESS + true + elseif result == LibGraphBLAS.GrB_NO_VALUE + false + else + @wraperror result + end + end + end +end + +# subassign fallback for Matrix <- Matrix, and Matrix <- Vector +""" + subassign!(C::GBMatrix, A::GBMatrix, I, J; kwargs...)::GBMatrix + +Assign a submatrix of `A` to `C`. Equivalent to [`assign!`](@ref) except that +`size(mask) == size(A)`, whereas `size(mask) == size(C)` in `assign!`. + +# Arguments +- `C::GBMatrix`: the matrix being subassigned to where `C[I,J] = A`. +- `A::GBMatrix`: the matrix being assigned to a submatrix of `C`. +- `I` and `J`: A colon, scalar, vector, or range indexing C. + +# Keywords +- `mask::Union{Nothing, GBMatrix} = nothing`: mask where + `size(M) == size(A)`. +- `accum::Union{Nothing, Function, AbstractBinaryOp} = 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 +- `GBMatrix`: The input matrix A. + +# Throws +- `GrB_DIMENSION_MISMATCH`: If `size(A) != (max(I), max(J))` or `size(A) != size(mask)`. +""" +function subassign!( + C::AbstractGBMatrix, A::GBArray, I, J; + mask = nothing, accum = nothing, desc = nothing +) + I, ni = idx(I) + J, nj = idx(J) + mask === nothing && (mask = C_NULL) + I = decrement!(I) + J = decrement!(J) + # we know A isn't adjoint/transpose on input + desc = _handledescriptor(desc; in1=A) + @wraperror LibGraphBLAS.GxB_Matrix_subassign(gbpointer(C), mask, getaccum(accum, eltype(C)), gbpointer(parent(A)), I, ni, J, nj, desc) + increment!(I) + increment!(J) + return A +end + +function subassign!(C::AbstractGBArray, x, I, J; + mask = nothing, accum = nothing, desc = nothing +) + I, ni = idx(I) + J, nj = idx(J) + I = decrement!(I) + J = decrement!(J) + desc = _handledescriptor(desc) + mask, accum = _handlenothings(mask, accum) + _subassign(C, x, I, ni, J, nj, mask, getaccum(accum, eltype(C)), desc) + increment!(I) + increment!(J) +end + +function subassign!(C::AbstractGBArray, x::AbstractArray, I, J; + mask = nothing, accum = nothing, desc = nothing) + as(GBMatrix, x) do array + subassign!(C, array, I, J; mask, accum, desc) + end +end + +""" + assign!(C::GBMatrix, A::GBMatrix, I, J; kwargs...)::GBMatrix + +Assign a submatrix of `A` to `C`. Equivalent to [`subassign!`](@ref) except that +`size(mask) == size(C)`, whereas `size(mask) == size(A) in `subassign!`. + +# Arguments +- `C::GBMatrix`: the matrix being subassigned to where `C[I,J] = A`. +- `A::GBMatrix`: the matrix being assigned to a submatrix of `C`. +- `I` and `J`: A colon, scalar, vector, or range indexing C. + +# Keywords +- `mask::Union{Nothing, GBMatrix} = nothing`: mask where + `size(M) == size(C)`. +- `accum::Union{Nothing, Function, AbstractBinaryOp} = 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 +- `GBMatrix`: The input matrix A. + +# Throws +- `GrB_DIMENSION_MISMATCH`: If `size(A) != (max(I), max(J))` or `size(C) != size(mask)`. +""" +function assign!( + C::AbstractGBMatrix, A::AbstractGBVector, I, J; + mask = nothing, accum = nothing, desc = nothing +) + I, ni = idx(I) + J, nj = idx(J) + mask === nothing && (mask = C_NULL) + 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) + increment!(I) + increment!(J) + return A +end + +function assign!(C::AbstractGBArray, x, I, J; + mask = nothing, accum = nothing, desc = nothing +) + I = decrement!(I) + J = decrement!(J) + desc = _handledescriptor(desc) + _assign(gbpointer(C), x, I, ni, J, nj, mask, getaccum(accum, eltype(C)), desc) + increment!(I) + increment!(J) +end + +# setindex! uses subassign rather than assign. +function Base.setindex!( + C::AbstractGBMatrix, A, ::Colon, J; + mask = nothing, accum = nothing, desc = nothing +) + subassign!(C, A, ALL, J; mask, accum, desc) +end +function Base.setindex!( + C::AbstractGBMatrix, A, I, ::Colon; + mask = nothing, accum = nothing, desc = nothing +) + subassign!(C, A, I, ALL; mask, accum, desc) +end +function Base.setindex!( + C::AbstractGBMatrix, A, ::Colon, ::Colon; + mask = nothing, accum = nothing, desc = nothing +) + subassign!(C, A, ALL, ALL; mask, accum, desc) +end + +function Base.setindex!( + C::AbstractGBMatrix, + A, + I::Union{Vector, UnitRange, StepRange, Number}, + J::Union{Vector, UnitRange, StepRange, Number}; + mask = nothing, + accum = nothing, + desc = nothing +) + subassign!(C, A, I, J; mask, accum, desc) +end + +#Help wanted: This isn't really centered for a lot of eltypes. +function Base.replace_in_print_matrix(A::AbstractGBMatrix, i::Integer, j::Integer, s::AbstractString) + Base.isstored(A, i, j) ? s : Base.replace_with_centered_mark(s) +end + +# AbstractGBVector functions: +############################# +function Base.size(v::AbstractGBVector) + nrows = Ref{LibGraphBLAS.GrB_Index}() + @wraperror LibGraphBLAS.GrB_Matrix_nrows(nrows, gbpointer(v)) + return (Int64(nrows[]),) +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) + return v +end + +function Base.resize!(v::AbstractGBVector, n) + @wraperror LibGraphBLAS.GrB_Matrix_resize(gbpointer(v), n, 1) + return v +end + +function LinearAlgebra.diag(A::AbstractGBMatrix{T}, k::Integer = 0; desc = nothing) where {T} + m, n = size(A) + if !(k in -m:n) + s = 0 + elseif k >= 0 + s = min(m, n - k) + else + s = min(m + k, n) + end + v = GBVector{T}(s; A.fill) + desc = _handledescriptor(desc; in1=A) + if A isa Transpose + k = -k + end + @wraperror LibGraphBLAS.GxB_Vector_diag(LibGraphBLAS.GrB_Vector(gbpointer(v)), gbpointer(parent(A)), k, desc) + return v +end + +# This does not conform to the normal definition with a lazy wrapper. +function LinearAlgebra.Diagonal(v::AbstractGBVector, k::Integer=0; desc = nothing) + s = size(v, 1) + C = GBMatrix{eltype(v)}(s, s; fill = v.fill) + desc = _handledescriptor(desc) + # Switch ptr to a Vector to trick GraphBLAS. + # This is allowed since GrB_Vector is a GrB_Matrix internally. + @wraperror LibGraphBLAS.GxB_Matrix_diag(C, Ptr{LibGraphBLAS.GrB_Vector}(gbpointer(v)), k, desc) + return C +end + +# Type dependent functions build, setindex, getindex, and findnz: +for T ∈ valid_vec + if T ∈ gxb_vec + prefix = :GxB + else + prefix = :GrB + end + # Build functions + func = Symbol(prefix, :_Matrix_build_, suffix(T)) + @eval begin + function build(v::AbstractGBVector{$T}, I::Vector{<:Integer}, X::Vector{$T}; combine = +) + nnz(v) == 0 || throw(OutputNotEmptyError("Cannot build vector with existing elements")) + 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) + decrement!(I) + @wraperror LibGraphBLAS.$func( + Ptr{LibGraphBLAS.GrB_Vector}(gbpointer(v)), + I, + # TODO, fix this ugliness by switching to the GBVector build internally. + zeros(LibGraphBLAS.GrB_Index, length(I)), + X, + length(X), + combine + ) + increment!(I) + end + end + # Setindex functions + func = Symbol(prefix, :_Matrix_setElement_, suffix(T)) + @eval begin + function Base.setindex!(v::AbstractGBVector{$T}, x, i::Integer) + x = convert($T, x) + return LibGraphBLAS.$func(gbpoitner(v), x, LibGraphBLAS.GrB_Index(decrement!(i)), 0) + end + end + # Getindex functions + func = Symbol(prefix, :_Matrix_extractElement_, suffix(T)) + @eval begin + function Base.getindex(v::GBVector{$T}, i::Integer) + x = Ref{$T}() + result = LibGraphBLAS.$func(x, v, LibGraphBLAS.GrB_Index(decrement!(i)), 0) + if result == LibGraphBLAS.GrB_SUCCESS + return x[] + elseif result == LibGraphBLAS.GrB_NO_VALUE + return v.fill + else + @wraperror result + end + end + end + # findnz functions + func = Symbol(prefix, :_Matrix_extractTuples_, suffix(T)) + @eval begin + function SparseArrays.findnz(v::AbstractGBVector{$T}) + nvals = Ref{LibGraphBLAS.GrB_Index}(nnz(v)) + I = Vector{LibGraphBLAS.GrB_Index}(undef, nvals[]) + X = Vector{$T}(undef, nvals[]) + wait(v) + @wraperror LibGraphBLAS.$func(I, C_NULL, X, nvals, gbpointer(v)) + nvals[] == length(I) == length(X) || throw(DimensionMismatch("length(I) != length(X)")) + return increment!(I), X + end + function SparseArrays.nonzeros(v::GBVector{$T}) + nvals = Ref{LibGraphBLAS.GrB_Index}(nnz(v)) + X = Vector{$T}(undef, nvals[]) + wait(v) + @wraperror LibGraphBLAS.$func(C_NULL, C_NULL, X, nvals, gbpointer(v)) + nvals[] == length(X) || throw(DimensionMismatch("")) + return X + end + function SparseArrays.nonzeroinds(v::GBVector{$T}) + nvals = Ref{LibGraphBLAS.GrB_Index}(nnz(v)) + I = Vector{LibGraphBLAS.GrB_Index}(undef, nvals[]) + wait(v) + @wraperror LibGraphBLAS.$func(I, C_NULL, C_NULL, nvals, gbpointer(v)) + nvals[] == length(I) || throw(DimensionMismatch("")) + return increment!(I) + end + end +end + +function build(v::GBVector{T}, I::Vector, x::T) where {T} + nnz(v) == 0 || throw(OutputNotEmptyError("Cannot build vector with existing elements")) + x = GBScalar(x) + decrement!(I) + @wraperror LibGraphBLAS.GxB_Matrix_build_Scalar( + v, + Vector{LibGraphBLAS.GrB_Index}(I), + zeros(LibGraphBLAS.GrB_Index, length(I)), + x, + length(I) + ) + increment!(I) + return v +end + +""" + subassign(w::GBVector, u::GBVector, I; kwargs...)::GBVector + +Assign a subvector of `w` to `u`. Return `u`. Equivalent to the matrix definition. +""" +function subassign!(w::AbstractGBVector{T, F}, u, I; mask = nothing, accum = nothing, desc = nothing) where {T, F} + return subassign!(GBMatrix{T, F}(w.p, w.fill), u, I, UInt64[1]; mask, accum, desc) +end + +""" + assign(w::GBVector, u::GBVector, I; kwargs...)::GBVector + +Assign a subvector of `w` to `u`. Return `u`. Equivalent to the matrix definition. +""" +function assign!(w::AbstractGBVector{T, F}, u, I; mask = nothing, accum = nothing, desc = nothing) where {T, F} + return assign!(GBMatrix{T, F}(w.p, w.fill), u, I, UInt64[1]; mask, accum, desc) +end + +function Base.setindex!( + u::AbstractGBVector, x, ::Colon; + mask = nothing, accum = nothing, desc = nothing +) + subassign!(u, x, ALL; mask, accum, desc) + return nothing +end +# silly overload to help a bit with broadcasting. +function Base.setindex!( + u::AbstractGBVector, x, I::Union{Vector, UnitRange, StepRange, Colon}, ::Colon; + mask = nothing, accum = nothing, desc = nothing +) + Base.setindex!(u, x, I; mask, accum, desc) +end +function Base.setindex!( + u::AbstractGBVector, x, I::Union{Vector, UnitRange, StepRange}; + mask = nothing, accum = nothing, desc = nothing +) + subassign!(u, x, I; mask, accum, desc) + return nothing +end \ No newline at end of file diff --git a/src/abstracts.jl b/src/abstracts.jl index 566a1b4d..f9d38a35 100644 --- a/src/abstracts.jl +++ b/src/abstracts.jl @@ -7,3 +7,8 @@ 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 + +abstract type AbstractGBMatrix{T, F} <: AbstractGBArray{T, 2, F} end +abstract type AbstractGBVector{T, F} <: AbstractGBArray{T, 1, F} end \ No newline at end of file diff --git a/src/asjulia.jl b/src/asjulia.jl index 26d5e289..3ae73735 100644 --- a/src/asjulia.jl +++ b/src/asjulia.jl @@ -1,6 +1,21 @@ -function as(f::Function, type::Type{<:Union{Matrix, Vector}}, A::GBVecOrMat{T}; dropzeros=false, freeunpacked=false, nomodstructure = false) where {T} - (type == Matrix && !(A isa GBMatrix)) && throw(ArgumentError("Cannot wrap $(typeof(A)) in a Matrix.")) - (type == Vector && !(A isa GBVector)) && throw(ArgumentError("Cannot wrap $(typeof(A)) in a Vector.")) +function as(f::Function, type::Type{<:Union{GBMatrix, GBVector}}, A::AbstractArray{T}) where {T} + if !(A isa DenseVecOrMat) + A = A isa AbstractVector ? collect(A) : Matrix(A) + end + if type <:AbstractGBMatrix + array = type{T}(size(A, 1), size(A, 2)) + else + array = type{T}(size(A, 1)) + end + _packdensematrix!(array, A) + _makeshallow!(array) + return f(array) +end + + +function as(f::Function, type::Type{<:Union{Matrix, Vector}}, A::AbstractGBArray{T}; dropzeros=false, freeunpacked=false, nomodstructure = false) where {T} + (type == Matrix && !(A isa AbstractMatrix)) && throw(ArgumentError("Cannot wrap $(typeof(A)) in a Matrix.")) + (type == Vector && !(A isa AbstractVector)) && throw(ArgumentError("Cannot wrap $(typeof(A)) in a Vector.")) if gbget(A, SPARSITY_STATUS) != Int64(GBDENSE) X = similar(A) if X isa GBVector @@ -31,7 +46,7 @@ function as(f::Function, type::Type{<:Union{Matrix, Vector}}, A::GBVecOrMat{T}; return result end -function as(f::Function, ::Type{SparseMatrixCSC}, A::GBMatrix{T}; freeunpacked=false) where {T} +function as(f::Function, ::Type{SparseMatrixCSC}, A::AbstractGBMatrix{T}; freeunpacked=false) where {T} colptr, colptrsize, rowidx, rowidxsize, values, valsize = _unpackcscmatrix!(A) array = SparseMatrixCSC{T, LibGraphBLAS.GrB_Index}(size(A, 1), size(A, 2), colptr, rowidx, values) result = try @@ -48,7 +63,7 @@ function as(f::Function, ::Type{SparseMatrixCSC}, A::GBMatrix{T}; freeunpacked=f return result end -function as(f::Function, ::Type{SparseVector}, A::GBVector{T}; freeunpacked=false) where {T} +function as(f::Function, ::Type{SparseVector}, A::AbstractGBVector{T}; freeunpacked=false) where {T} colptr, colptrsize, rowidx, rowidxsize, values, valsize = _unpackcscmatrix!(A) vector = SparseVector{T, LibGraphBLAS.GrB_Index}(size(A, 1), rowidx, values) result = try @@ -66,41 +81,41 @@ function as(f::Function, ::Type{SparseVector}, A::GBVector{T}; freeunpacked=fals end -function Base.Matrix(A::GBMatrix) +function Base.Matrix(A::AbstractGBMatrix) # we use nomodstructure here to avoid the pitfall of densifying A. return as(Matrix, A; nomodstructure=true) do arr return copy(arr) end end -function Matrix!(A::GBMatrix) +function Matrix!(A::AbstractGBMatrix) # we use nomodstructure here to avoid the pitfall of densifying A. return as(Matrix, A; freeunpacked=true) do arr return copy(arr) end end -function Base.Vector(v::GBVector) +function Base.Vector(v::AbstractGBVector) # we use nomodstructure here to avoid the pitfall of densifying A. return as(Vector, v; nomodstructure=true) do vec return copy(vec) end end -function Vector!(v::GBVector) +function Vector!(v::AbstractGBVector) # we use nomodstructure here to avoid the pitfall of densifying A. return as(Vector, v; freeunpacked=true) do vec return copy(vec) end end -function SparseArrays.SparseMatrixCSC(A::GBMatrix) +function SparseArrays.SparseMatrixCSC(A::AbstractGBMatrix) return as(SparseMatrixCSC, A) do arr return copy(arr) end end -function SparseArrays.SparseVector(v::GBVector) +function SparseArrays.SparseVector(v::AbstractGBVector) return as(SparseVector, v) do arr return copy(arr) end diff --git a/src/constants.jl b/src/constants.jl index 59555ea0..e8318ecb 100644 --- a/src/constants.jl +++ b/src/constants.jl @@ -1,7 +1,12 @@ -const GBVecOrMat{T} = Union{GBVector{T}, GBMatrix{T}} -const GBMatOrTranspose{T} = Union{GBMatrix{T}, Transpose{T, GBMatrix{T}}} -const GBVecOrTranspose{T} = Union{GBVector{T}, Transpose{T, GBVector{T}}} -const GBArray{T} = Union{GBVecOrTranspose{T}, GBMatOrTranspose{T}} +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 ptrtogbtype = IdDict{Ptr, GBType}() const GrBOp = Union{ diff --git a/src/import.jl b/src/import.jl index 7110b065..769f184f 100644 --- a/src/import.jl +++ b/src/import.jl @@ -30,7 +30,7 @@ function _importcscmat( jumbled, desc ) - return A[] + return A end function _importcscmat( @@ -43,20 +43,19 @@ function _importcscmat( desc = nothing, iso = false ) where {U, T} - # This section comes after some chatting with Keno Fisher. - # Cannot directly pass Julia arrays to GraphBLAS, it expects malloc'd arrays. - # Instead we'll malloc some memory for each of the three vectors, and unsafe_copyto! - # into them. colsize = LibGraphBLAS.GrB_Index(sizeof(colptr)) #Size of colptr vector rowsize = LibGraphBLAS.GrB_Index(sizeof(rowindices)) #Size of rowindex vector valsize = LibGraphBLAS.GrB_Index(sizeof(values)) #Size of nzval vector col = ccall(:jl_malloc, Ptr{LibGraphBLAS.GrB_Index}, (UInt, ), colsize) - unsafe_copyto!(col, Ptr{UInt64}(pointer(colptr .- 1)), length(colptr)) + unsafe_copyto!(col, Ptr{UInt64}(pointer(decrement!(colptr))), length(colptr)) row = ccall(:jl_malloc, Ptr{LibGraphBLAS.GrB_Index}, (UInt, ), rowsize) - unsafe_copyto!(row, Ptr{UInt64}(pointer(rowindices .- 1)), length(rowindices)) + unsafe_copyto!(row, Ptr{UInt64}(pointer(decrement!(rowindices))), length(rowindices)) val = ccall(:jl_malloc, Ptr{T}, (UInt, ), valsize) unsafe_copyto!(val, pointer(values), length(values)) - return _importcscmat(m, n, col, colsize, row, rowsize, val, valsize; jumbled, desc, iso) + x = _importcscmat(m, n, col, colsize, row, rowsize, val, valsize; jumbled, desc, iso) + increment!(colptr) + increment!(rowindices) + return x end """ @@ -67,10 +66,11 @@ Create a GBMatrix from a SparseArrays.SparseMatrixCSC `S`. Note, that unlike other methods of construction, the resulting matrix will be held by column. Use `gbset(A, :format, :byrow)` to switch to row orientation. """ -function GBMatrix(S::SparseMatrixCSC) - return GBMatrix{eltype(S)}(_importcscmat(S.m, S.n, S.colptr, S.rowval, S.nzval)) +function GBMatrix(S::SparseMatrixCSC{T}; fill::F = nothing) where {T, F} + return GBMatrix{T, F}(_importcscmat(S.m, S.n, S.colptr, S.rowval, S.nzval), fill) end +# TODO: should be able to do better here. function GBMatrix(v::SparseVector) S = SparseMatrixCSC(v) return GBMatrix(S) @@ -81,44 +81,8 @@ end Create a GBVector from SparseArrays sparse vector `v`. """ -function GBVector(v::SparseVector) - return GBVector{eltype(v)}(_importcscmat(v.n, 1, [1, length(v.nzind) + 1], v.nzind, v.nzval)); -end - -function _importcscvec( - n::Integer, vi::Ptr{UInt64}, vi_size, vx::Ptr{T}, vx_size, nnz; - jumbled::Bool = false, desc = nothing, iso = false -) where {T} - v = Ref{LibGraphBLAS.GrB_Vector}() - n = LibGraphBLAS.GrB_Index(n) - desc = _handledescriptor(desc) - @wraperror LibGraphBLAS.GxB_Vector_import_CSC( - v, - gbtype(T), - n, - Ref{Ptr{LibGraphBLAS.GrB_Index}}(vi), - Ref{Ptr{Cvoid}}(vx), - vi_size, - vx_size, - iso, - nnz, - jumbled, - desc - ) - return GBVector{T}(v[]) -end - -function _importcscvec( - n::Integer, vi::Vector{U}, vx::Vector{T}, nnz; - jumbled::Bool = false, desc = nothing, iso = false -) where {U,T} - vi_size = LibGraphBLAS.GrB_Index(sizeof(vi)) - vx_size = LibGraphBLAS.GrB_Index(sizeof(vx)) - indices = ccall(:jl_malloc, Ptr{LibGraphBLAS.GrB_Index}, (UInt, ), vi_size) - unsafe_copyto!(indices, Ptr{UInt64}(pointer(vi .- 1)), length(vi)) - values = ccall(:jl_malloc, Ptr{T}, (UInt, ), vx_size) - unsafe_copyto!(values, pointer(vx), length(vx)) - return _importcscvec(n, indices, vi_size, values, vx_size, nnz; jumbled, desc, iso) +function GBVector(v::SparseVector{T}; fill::F = nothing) where {T, F} + return GBVector{T, F}(_importcscmat(v.n, 1, [1, length(v.nzind) + 1], v.nzind, v.nzval), fill) end function _importcsrmat( @@ -153,7 +117,7 @@ function _importcsrmat( jumbled, desc ) - return GBMatrix{T}(A[]) + return A end function _importcsrmat( @@ -204,7 +168,7 @@ function _importdensematrix( iso, desc ) - return C[] + return C end function _importdensematrix( @@ -224,56 +188,24 @@ end Create a GBMatrix from a Julia dense matrix. """ -function GBMatrix(M::Union{AbstractVector, AbstractMatrix}) - #if M isa AbstractVector && !(M isa Vector) - # M = Vector(M) - #end - #if M isa AbstractMatrix && !(M isa Matrix) - # M = Matrix(M) - #end - return GBMatrix{eltype(M)}(_importdensematrix(size(M, 1), size(M, 2), M)) +function GBMatrix(M::Union{AbstractVector{T}, AbstractMatrix{T}}; fill::F = nothing) where {T, F} + if M isa AbstractVector && !(M isa Vector) + M = collect(M) + end + if M isa AbstractMatrix && !(M isa Matrix) + M = Matrix(M) + end + return GBMatrix{T, F}(_importdensematrix(size(M, 1), size(M, 2), M), fill) end """ - GBVector(v::SparseVector) + GBVector(v::AbstractVector) Create a GBVector from a Julia dense vector. """ -function GBVector(v::AbstractVector) +function GBVector(v::AbstractVector{T}; fill::F = nothing) where {T, F} if !(v isa Vector) - v = Vector(v) + v = collect(v) end - return GBVector{eltype(v)}(_importdensematrix(size(v, 1), 1, v)) -end - -function _importdensevec( - n::Integer, v::Ptr{T}, vsize; - desc = nothing, iso = false -) where {T} - w = Ref{LibGraphBLAS.GrB_Vector}() - n = LibGraphBLAS.GrB_Index(n) - desc = _handledescriptor(desc) - @wraperror LibGraphBLAS.GxB_Vector_import_Full( - w, - gbtype(T), - n, - Ref{Ptr{Cvoid}}(v), - vsize, - iso, - desc - ) - wout = GBVector{T}(w[]) - return wout -end - -function _importdensevec( - n::Integer, v::Vector{T}; - desc = nothing, iso = false -) where {T} - n = LibGraphBLAS.GrB_Index(n) - vsize = LibGraphBLAS.GrB_Index(sizeof(v)) - # We have to do this instead of Libc.malloc because GraphBLAS will use :jl_free, not Libc.free - vx = ccall(:jl_malloc, Ptr{T}, (UInt, ), vsize) - unsafe_copyto!(vx, pointer(v), length(v)) - return _importdensevec(n, vx, vsize; desc, iso) -end + return GBVector{T, F}(_importdensematrix(size(v, 1), 1, v), fill) +end \ No newline at end of file diff --git a/src/indexutils.jl b/src/indexutils.jl index 53bddbc3..e261284e 100644 --- a/src/indexutils.jl +++ b/src/indexutils.jl @@ -6,6 +6,7 @@ proper format for GraphBLAS indexing. Should *not* be used for functions that ta scalar index like [`extractElement`]. """ function idx(I) + # TODO. Do better here, and minimize manual idx management in rest of library. if I == ALL return I, 0 #ni doesn't matter if I=ALL elseif I isa UnitRange @@ -19,7 +20,7 @@ function idx(I) return [I.start, I.stop, -I.step + 1], LibGraphBLAS.GxB_BACKWARDS end elseif I isa Vector - return Vector{LibGraphBLAS.GrB_Index}(I), length(I) #Assume ni = length(I) otherwise + return reinterpret(LibGraphBLAS.GrB_Index, I), length(I) #Assume ni = length(I) otherwise elseif I isa Integer return [UInt64(I)], 1 elseif I isa CartesianIndex{1} diff --git a/src/matrix.jl b/src/matrix.jl index 2a0805b2..42b1b019 100644 --- a/src/matrix.jl +++ b/src/matrix.jl @@ -1,20 +1,22 @@ # Constructors: ############### """ - GBMatrix{T}(nrows = LibGraphBLAS.GxB_INDEX_MAX, ncols = LibGraphBLAS.GxB_INDEX_MAX) + GBMatrix{T}(nrows, ncols; fill = nothing) Create a GBMatrix of the specified size, defaulting to the maximum on each dimension, 2^60. """ -function GBMatrix{T}(nrows::Integer = LibGraphBLAS.GxB_INDEX_MAX, ncols::Integer = LibGraphBLAS.GxB_INDEX_MAX) where {T} +function GBMatrix{T}(nrows::Integer, ncols::Integer; fill::F = nothing) where {T, F} m = Ref{LibGraphBLAS.GrB_Matrix}() @wraperror LibGraphBLAS.GrB_Matrix_new(m, gbtype(T), nrows, ncols) - return GBMatrix{T}(m[]) + return GBMatrix{T, F}(finalizer(m) do ref + @wraperror LibGraphBLAS.GrB_Matrix_free(ref) + end, fill) end -GBMatrix{T}(dims::Dims{2}) where {T} = GBMatrix{T}(dims...) -GBMatrix{T}(dims::Tuple{<:Integer}) where {T} = GBMatrix{T}(dims...) -GBMatrix{T}(size::Tuple{Base.OneTo, Base.OneTo}) where {T} = - GBMatrix{T}(size[1].stop, size[2].stop) +GBMatrix{T}(dims::Dims{2}; fill = nothing) where {T} = GBMatrix{T}(dims...; fill) +GBMatrix{T}(dims::Tuple{<:Integer}; fill = nothing) where {T} = GBMatrix{T}(dims...; fill) +GBMatrix{T}(size::Tuple{Base.OneTo, Base.OneTo}; fill = nothing) where {T} = + GBMatrix{T}(size[1].stop, size[2].stop; fill) """ GBMatrix(I, J, X; combine = +, nrows = maximum(I), ncols = maximum(J)) @@ -24,9 +26,12 @@ to `|` for booleans and `+` for nonbooleans. """ function GBMatrix( I::AbstractVector, J::AbstractVector, X::AbstractVector{T}; - combine = +, nrows = maximum(I), ncols = maximum(J) + combine = +, nrows = maximum(I), ncols = maximum(J), fill = nothing ) where {T} - A = GBMatrix{T}(nrows, ncols) + I isa Vector || (I = collect(I)) + J isa Vector || (J = collect(J)) + X isa Vector || (X = collect(X)) + A = GBMatrix{T}(nrows, ncols; fill) build(A, I, J, X; combine) return A end @@ -40,21 +45,22 @@ The resulting matrix is "iso-valued" such that it only stores `x` once rather th each index. """ function GBMatrix(I::AbstractVector, J::AbstractVector, x::T; - nrows = maximum(I), ncols = maximum(J)) where {T} - A = GBMatrix{T}(nrows, ncols) + nrows = maximum(I), ncols = maximum(J), fill = nothing) where {T} + A = GBMatrix{T}(nrows, ncols; fill) build(A, I, J, x) return A end -function GBMatrix(dims::Dims{2}, x::T) where {T} - A = GBMatrix{T}(dims) +function GBMatrix(dims::Dims{2}, x::T; fill = nothing) where {T} + A = GBMatrix{T}(dims; fill) A[:, :] = x return A end -GBMatrix(nrows, ncols, x::T) where {T} = GBMatrix((nrows, ncols), x) +GBMatrix(nrows, ncols, x::T; fill::F = nothing) where {T, F} = GBMatrix((nrows, ncols), x; fill) +# TODO, switch to pointer fudging. function GBMatrix(v::GBVector) A = GBMatrix{eltype(v)}(size(v, 1), size(v, 2)) nz = findnz(v) @@ -64,85 +70,60 @@ function GBMatrix(v::GBVector) return A end -function build(A::GBMatrix{T}, I::AbstractVector, J::AbstractVector, x::T) where {T} - nnz(A) == 0 || throw(OutputNotEmptyError("Cannot build matrix with existing elements")) - length(I) == length(J) || DimensionMismatch("I, J and X must have the same length") - x = GBScalar(x) - - @wraperror LibGraphBLAS.GxB_Matrix_build_Scalar( - A, - Vector{LibGraphBLAS.GrB_Index}(decrement!(I)), - Vector{LibGraphBLAS.GrB_Index}(decrement!(J)), - x, - length(I) - ) - increment!(I) - increment!(J) - return A -end - - # Some Base and basic SparseArrays/LinearAlgebra functions: ########################################################### -Base.unsafe_convert(::Type{LibGraphBLAS.GrB_Matrix}, A::GBMatrix) = A.p +Base.unsafe_convert(::Type{LibGraphBLAS.GrB_Matrix}, A::GBMatrix) = A.p[] -function Base.copy(A::GBMatrix{T}) where {T} +function Base.copy(A::GBMatrix{T, F}) where {T, F} C = Ref{LibGraphBLAS.GrB_Matrix}() - LibGraphBLAS.GrB_Matrix_dup(C, A) - return GBMatrix{T}(C[]) + LibGraphBLAS.GrB_Matrix_dup(C, gbpointer(A)) + return GBMatrix{T, F}(C, A.fill) end -""" - clear!(v::GBVector) - clear!(A::GBMatrix) -Clear all the entries from the GBArray. -Does not modify the type or dimensions. -""" -clear!(A::GBArray) = @wraperror LibGraphBLAS.GrB_Matrix_clear(parent(A)); return nothing - -function Base.size(A::GBMatrix) - nrows = Ref{LibGraphBLAS.GrB_Index}() - ncols = Ref{LibGraphBLAS.GrB_Index}() - @wraperror LibGraphBLAS.GrB_Matrix_nrows(nrows, A) - @wraperror LibGraphBLAS.GrB_Matrix_ncols(ncols, A) - return (Int64(nrows[]), Int64(ncols[])) +# because of the fill kwarg we have to redo a lot of the Base.similar dispatch stack. +function Base.similar( + A::GBMatrixOrTranspose{T}, ::Type{TNew} = T, + dims::Tuple{Int64, Vararg{Int64, N}} = size(A); fill = parent(A).fill +) where {T, TNew, N} + if dims isa Dims{1} + return GBVector{TNew}(dims...; fill) + else + return GBMatrix{TNew}(dims...; fill) + end end -function SparseArrays.nnz(A::GBArray) - nvals = Ref{LibGraphBLAS.GrB_Index}() - @wraperror LibGraphBLAS.GrB_Matrix_nvals(nvals, parent(A)) - return Int64(nvals[]) +function Base.similar(A::GBMatrixOrTranspose{T}, dims::Tuple; fill = parent(A).fill) where {T} + return similar(A, T, dims; fill) end -Base.eltype(::Type{GBMatrix{T}}) where{T} = T - function Base.similar( - ::GBMatrix{T}, ::Type{TNew}, - dims::Union{Dims{1}, Dims{2}} + A::GBMatrixOrTranspose{T}, ::Type{TNew}, + dims::Integer; fill = parent(A).fill ) where {T, TNew} - return GBMatrix{TNew}(dims...) + return similar(A, TNew, (dims,); fill) end -function Base.deleteat!(A::GBMatrix, i, j) - @wraperror LibGraphBLAS.GrB_Matrix_removeElement(A, decrement!(i), decrement!(j)) - return A +function Base.similar( + A::GBMatrixOrTranspose{T}, ::Type{TNew}, + dim1::Integer, dim2::Integer; fill = parent(A).fill +) where {T, TNew} + return similar(A, TNew, (dim1, dim2); fill) end -function Base.resize!(A::GBMatrix, nrows_new, ncols_new) - @wraperror LibGraphBLAS.GrB_Matrix_resize(A, nrows_new, ncols_new) - return A +function Base.similar( + A::GBMatrixOrTranspose{T}, + dims::Integer; fill = parent(A).fill +) where {T} + return similar(A, (dims,); fill) end -# This does not conform to the normal definition with a lazy wrapper. -function LinearAlgebra.Diagonal(v::GBVector, k::Integer=0; desc = nothing) - s = size(v, 1) - C = GBMatrix{eltype(v)}(s, s) - desc = _handledescriptor(desc) - # Switch ptr to a Vector to trick GraphBLAS. - # This is allowed since GrB_Vector is a GrB_Matrix internally. - @wraperror LibGraphBLAS.GxB_Matrix_diag(C, Ptr{LibGraphBLAS.GrB_Vector}(v.p), k, desc) - return C + +function Base.similar( + A::GBMatrixOrTranspose{T}, + dim1::Integer, dim2::Integer; fill = parent(A).fill +) where {T} + return similar(A, (dim1, dim2); fill) end # TODO: FIXME @@ -150,103 +131,6 @@ end # return Diagonal(v, k; desc) # end -# Type dependent functions build, setindex, getindex, and findnz: -for T ∈ valid_vec - if T ∈ gxb_vec - prefix = :GxB - else - prefix = :GrB - end - # Build functions - func = Symbol(prefix, :_Matrix_build_, suffix(T)) - @eval begin - function build(A::GBMatrix{$T}, I::AbstractVector, J::AbstractVector, X::Vector{$T}; - combine = + - ) - combine = BinaryOp(combine)($T) - if !(I isa Vector) - I = Vector(I) - end - if !(J isa Vector) - J = Vector(J) - end - nnz(A) == 0 || throw(OutputNotEmptyError("Cannot build matrix with existing elements")) - length(X) == length(I) == length(J) || - DimensionMismatch("I, J and X must have the same length") - decrement!(I) - decrement!(J) - @wraperror LibGraphBLAS.$func( - A, - Vector{LibGraphBLAS.GrB_Index}(I), - Vector{LibGraphBLAS.GrB_Index}(J), - X, - length(X), - combine - ) - increment!(I) - increment!(J) - end - end - # Setindex functions - func = Symbol(prefix, :_Matrix_setElement_, suffix(T)) - @eval begin - function Base.setindex!(A::GBMatrix{$T}, x, i::Integer, j::Integer) - x = convert($T, x) - @wraperror LibGraphBLAS.$func(A, x, LibGraphBLAS.GrB_Index(decrement!(i)), LibGraphBLAS.GrB_Index(decrement!(j))) - return x - end - end - # Getindex functions - func = Symbol(prefix, :_Matrix_extractElement_, suffix(T)) - @eval begin - function Base.getindex(A::GBMatrix{$T}, i::Int, j::Int) - x = Ref{$T}() - result = LibGraphBLAS.$func(x, A, decrement!(i), decrement!(j)) - if result == LibGraphBLAS.GrB_SUCCESS - return x[] - elseif result == LibGraphBLAS.GrB_NO_VALUE - return nothing - else - @wraperror result - end - end - # Fix ambiguity - function Base.getindex(A::Transpose{$T, GBMatrix{$T}}, i::Int, j::Int) - return getindex(parent(A), j, i) - end - end - # findnz functions - func = Symbol(prefix, :_Matrix_extractTuples_, suffix(T)) - @eval begin - function SparseArrays.findnz(A::GBMatrix{$T}) - nvals = Ref{LibGraphBLAS.GrB_Index}(nnz(A)) - I = Vector{LibGraphBLAS.GrB_Index}(undef, nvals[]) - J = Vector{LibGraphBLAS.GrB_Index}(undef, nvals[]) - X = Vector{$T}(undef, nvals[]) - @wraperror LibGraphBLAS.$func(I, J, X, nvals, A) - nvals[] == length(I) == length(J) == length(X) || throw(DimensionMismatch("length(I) != length(X)")) - return increment!(I), increment!(J), X - end - function SparseArrays.nonzeros(A::GBMatrix{$T}) - nvals = Ref{LibGraphBLAS.GrB_Index}(nnz(A)) - X = Vector{$T}(undef, nvals[]) - @wraperror LibGraphBLAS.$func(C_NULL, C_NULL, X, nvals, A) - nvals[] == length(X) || throw(DimensionMismatch("")) - return X - end - function SparseArrays.nonzeroinds(A::GBMatrix{$T}) - nvals = Ref{LibGraphBLAS.GrB_Index}(nnz(A)) - I = Vector{LibGraphBLAS.GrB_Index}(undef, nvals[]) - J = Vector{LibGraphBLAS.GrB_Index}(undef, nvals[]) - wait(A) - @wraperror LibGraphBLAS.$func(I, J, C_NULL, nvals, A) - nvals[] == length(I) == length(J) || throw(DimensionMismatch("")) - return increment!(I), increment!(J) - end - end -end - - function Base.show(io::IO, ::MIME"text/plain", A::GBMatrix) gxbprint(io, A) end @@ -279,181 +163,3 @@ function Base.getindex(A::GBMatOrTranspose, v::AbstractVector) throw("Not implemented") end -""" - subassign!(C::GBMatrix, A::GBMatrix, I, J; kwargs...)::GBMatrix - -Assign a submatrix of `A` to `C`. Equivalent to [`assign!`](@ref) except that -`size(mask) == size(A)`, whereas `size(mask) == size(C)` in `assign!`. - -# Arguments -- `C::GBMatrix`: the matrix being subassigned to where `C[I,J] = A`. -- `A::GBMatrix`: the matrix being assigned to a submatrix of `C`. -- `I` and `J`: A colon, scalar, vector, or range indexing C. - -# Keywords -- `mask::Union{Nothing, GBMatrix} = nothing`: mask where - `size(M) == size(A)`. -- `accum::Union{Nothing, Function, AbstractBinaryOp} = 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 -- `GBMatrix`: The input matrix A. - -# Throws -- `GrB_DIMENSION_MISMATCH`: If `size(A) != (max(I), max(J))` or `size(A) != size(mask)`. -""" - -for T ∈ valid_vec - func = Symbol(:GxB_Matrix_subassign_, suffix(T)) - @eval begin - function _subassign(C::GBMatrix{$T}, x, I, ni, J, nj, mask, accum, desc) - @wraperror LibGraphBLAS.$func(C, mask, accum, x, I, ni, J, nj, desc) - return x - end - end - if T ∈ gxb_vec - prefix = :GxB - else - prefix = :GrB - end - func = Symbol(prefix, :_Matrix_assign_, suffix(T)) - @eval begin - function _assign(C::GBMatrix{$T}, x, I, ni, J, nj, mask, accum, desc) - @wraperror LibGraphBLAS.$func(C, mask, accum, x, I, ni, J, nj, desc) - return x - end - end -end -function subassign!( - C::GBMatrix, A, I, J; - mask = nothing, accum = nothing, desc = nothing -) - A1 = A - I, ni = idx(I) - J, nj = idx(J) - if A isa GBArray - elseif A isa AbstractVector - A = GBVector(A) - elseif A isa AbstractMatrix - A = GBMatrix(A) - end - mask === nothing && (mask = C_NULL) - I = decrement!(I) - J = decrement!(J) - if A isa GBArray - desc = _handledescriptor(desc; in1 = A) - @wraperror LibGraphBLAS.GxB_Matrix_subassign(C, mask, getaccum(accum, eltype(C)), parent(A), I, ni, J, nj, desc) - else - desc = _handledescriptor(desc) - _subassign(C, A, I, ni, J, nj, mask, getaccum(accum, eltype(C)), desc) - end - increment!(I) - increment!(J) - return A1 -end - -""" - assign!(C::GBMatrix, A::GBMatrix, I, J; kwargs...)::GBMatrix - -Assign a submatrix of `A` to `C`. Equivalent to [`subassign!`](@ref) except that -`size(mask) == size(C)`, whereas `size(mask) == size(A) in `subassign!`. - -# Arguments -- `C::GBMatrix`: the matrix being subassigned to where `C[I,J] = A`. -- `A::GBMatrix`: the matrix being assigned to a submatrix of `C`. -- `I` and `J`: A colon, scalar, vector, or range indexing C. - -# Keywords -- `mask::Union{Nothing, GBMatrix} = nothing`: mask where - `size(M) == size(C)`. -- `accum::Union{Nothing, Function, AbstractBinaryOp} = 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 -- `GBMatrix`: The input matrix A. - -# Throws -- `GrB_DIMENSION_MISMATCH`: If `size(A) != (max(I), max(J))` or `size(C) != size(mask)`. -""" -function assign!( - C::GBMatrix, A, I, J; - mask = nothing, accum = nothing, desc = nothing -) - I, ni = idx(I) - J, nj = idx(J) - if A isa GBArray - elseif A isa AbstractVector - A = GBVector(A) - elseif A isa AbstractMatrix - A = GBMatrix(A) - end - I = decrement!(I) - J = decrement!(J) - mask === nothing && (mask = C_NULL) - if A isa GBArray - desc = _handledescriptor(desc; in1 = A) - @wraperror LibGraphBLAS.GrB_Matrix_assign(C, mask, getaccum(accum, eltype(C)), parent(A), I, ni, J, nj, desc) - else - desc = _handledescriptor(desc) - _assign(C, A, I, ni, J, nj, mask, getaccum(accum, eltype(C)), desc) - end - increment!(I) - increment!(J) - return A -end - -# setindex! uses subassign rather than assign. -function Base.setindex!( - C::GBMatrix, A, ::Colon, J; - mask = nothing, accum = nothing, desc = nothing -) - subassign!(C, A, ALL, J; mask, accum, desc) -end -function Base.setindex!( - C::GBMatrix, A, I, ::Colon; - mask = nothing, accum = nothing, desc = nothing -) - subassign!(C, A, I, ALL; mask, accum, desc) -end -function Base.setindex!( - C::GBMatrix, A, ::Colon, ::Colon; - mask = nothing, accum = nothing, desc = nothing -) - subassign!(C, A, ALL, ALL; mask, accum, desc) -end - -function Base.setindex!( - C::GBMatrix, - A, - I::Union{Vector, UnitRange, StepRange, Number}, - J::Union{Vector, UnitRange, StepRange, Number}; - mask = nothing, - accum = nothing, - desc = nothing -) - subassign!(C, A, I, J; mask, accum, desc) -end - -function Base.setindex!( - ::GBMatrix, A, ::AbstractVector; - mask = nothing, accum = nothing, desc = nothing -) - throw("Not implemented") -end - -#Printing fixes: -function Base.isstored(A::GBArray, i::Integer, j::Integer) - @boundscheck checkbounds(A, i, j) - if A[i, j] === nothing - return false - else - return true - end -end - -#Help wanted: This isn't really centered for a lot of eltypes. -function Base.replace_in_print_matrix(A::GBArray, i::Integer, j::Integer, s::AbstractString) - Base.isstored(A, i, j) ? s : Base.replace_with_centered_mark(s) -end diff --git a/src/operations/concat.jl b/src/operations/concat.jl index 86775ba4..f7f35597 100644 --- a/src/operations/concat.jl +++ b/src/operations/concat.jl @@ -1,6 +1,6 @@ function cat!(C::GBArray, tiles::AbstractArray{T}) where {T<:GBArray} tiles = permutedims(tiles) - @wraperror LibGraphBLAS.GxB_Matrix_concat(C, tiles, size(tiles,2), size(tiles,1), C_NULL) + @wraperror LibGraphBLAS.GxB_Matrix_concat(gbpointer(C), gbpointer.(tiles), size(tiles,2), size(tiles,1), C_NULL) return C end @@ -10,7 +10,15 @@ end Create a new array formed from the contents of `tiles` in the sense of a block matrix This doesn't exactly match the Julia notion of `cat`. """ -function Base.cat(tiles::VecOrMat{T}) where {T<:GBArray} +function Base.cat(tiles::VecOrMat{<:Union{AbstractGBMatrix{T, F}, AbstractGBVector{T, F}}}) where {T, F} + fills = getproperty.(tiles, :fill) + if F <: Union{Nothing, Missing} + fill = F() + elseif all(y->y==fills[1], fills) + fill = fills[1] + else + fill = zero(F) + end ncols = sum(size.(tiles[1,:], 2)) nrows = sum(size.(tiles[:, 1], 1)) types = eltype.(tiles) @@ -18,11 +26,9 @@ function Base.cat(tiles::VecOrMat{T}) where {T<:GBArray} for type ∈ types[2:end] t = promote_type(t, type) end - if tiles isa AbstractArray{<:GBVector} && ncols == 1 - C = GBVector{t}(nrows) - else - C = GBMatrix{t}(nrows,ncols) - end + sz = (tiles isa AbstractArray && ncols == 1) ? (nrows,) : (nrows, ncols) + + C = similar(tiles[1], t, sz; fill) return cat!(C, tiles) end diff --git a/src/operations/ewise.jl b/src/operations/ewise.jl index 7af2adae..2dc38b4e 100644 --- a/src/operations/ewise.jl +++ b/src/operations/ewise.jl @@ -35,19 +35,12 @@ function emul!( size(C) == size(A) == size(B) || throw(DimensionMismatch()) op = BinaryOp(op)(eltype(A), eltype(B)) accum = getaccum(accum, eltype(C)) - if op isa TypedSemiring - @wraperror LibGraphBLAS.GrB_Matrix_eWiseMult_Semiring(C, mask, accum, op, parent(A), parent(B), desc) - return C - elseif op isa TypedMonoid - @wraperror LibGraphBLAS.GrB_Matrix_eWiseMult_Monoid(C, mask, accum, op, parent(A), parent(B), desc) - return C - elseif op isa TypedBinaryOperator - @wraperror LibGraphBLAS.GrB_Matrix_eWiseMult_BinaryOp(C, mask, accum, op, parent(A), parent(B), desc) + if op isa TypedBinaryOperator + @wraperror LibGraphBLAS.GrB_Matrix_eWiseMult_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 operator.")) end - return C end """ @@ -83,11 +76,7 @@ function emul( desc = nothing ) t = inferbinarytype(eltype(A), eltype(B), op) - if A isa GBVector && B isa GBVector - C = GBVector{t}(size(A)) - else - C = GBMatrix{t}(size(A)) - end + C = similar(A, t, size(A); fill=_promotefill(parent(A).fill, parent(B).fill)) return emul!(C, A, B, op; mask, accum, desc) end @@ -130,19 +119,12 @@ function eadd!( size(C) == size(A) == size(B) || throw(DimensionMismatch()) op = BinaryOp(op)(eltype(A), eltype(B)) accum = getaccum(accum, eltype(C)) - if op isa TypedSemiring - @wraperror LibGraphBLAS.GrB_Matrix_eWiseAdd_Semiring(C, mask, accum, op, parent(A), parent(B), desc) - return C - elseif op isa TypedMonoid - @wraperror LibGraphBLAS.GrB_Matrix_eWiseAdd_Monoid(C, mask, accum, op, parent(A), parent(B), desc) - return C - elseif op isa TypedBinaryOperator - @wraperror LibGraphBLAS.GrB_Matrix_eWiseAdd_BinaryOp(C, mask, accum, op, parent(A), parent(B), desc) + 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.")) end - return C end """ @@ -177,11 +159,7 @@ function eadd( desc = nothing ) t = inferbinarytype(eltype(A), eltype(B), op) - if A isa GBVector && B isa GBVector - C = GBVector{t}(size(A)) - else - C = GBMatrix{t}(size(A)) - end + C = similar(A, t, size(A); fill=_promotefill(parent(A).fill, parent(B).fill)) return eadd!(C, A, B, op; mask, accum, desc) end @@ -226,12 +204,11 @@ function eunion!( op = BinaryOp(op)(eltype(A), eltype(B)) accum = getaccum(accum, eltype(C)) if op isa TypedBinaryOperator - @wraperror LibGraphBLAS.GxB_Matrix_eWiseUnion(C, mask, accum, op, parent(A), GBScalar(α), parent(B), GBScalar(β), desc) + @wraperror LibGraphBLAS.GxB_Matrix_eWiseUnion(gbpointer(C), mask, accum, op, gbpointer(parent(A)), GBScalar(α), gbpointer(parent(B)), GBScalar(β), desc) return C else throw(ArgumentError("$op is not a valid binary op.")) end - return C end """ @@ -266,39 +243,10 @@ function eunion( desc = nothing ) where {T, U} t = inferbinarytype(eltype(A), eltype(B), op) - if A isa GBVector && B isa GBVector - C = GBVector{t}(size(A)) - else - C = GBMatrix{t}(size(A)) - end + C = similar(A, t, size(A); fill=_promotefill(parent(A).fill, parent(B).fill)) return eunion!(C, A, α, B, β, op; mask, accum, desc) end - -# function emul!(C, A, B, op::Function; mask = nothing, accum = nothing, desc = nothing) -# emul!(C, A, B, BinaryOp(op); mask, accum, desc) -# end -# -# function emul(A, B, op::Function; mask = nothing, accum = nothing, desc = nothing) -# emul(A, B, BinaryOp(op); mask, accum, desc) -# end -# -# function eadd!(C, A, B, op::Function; mask = nothing, accum = nothing, desc = nothing) -# eadd!(C, A, B, BinaryOp(op); mask, accum, desc) -# end -# -# function eadd(A, B, op::Function; mask = nothing, accum = nothing, desc = nothing) -# eadd(A, B, BinaryOp(op); mask, accum, desc) -# end -# -# function eunion!(C, A, α, B, β, op::Function; mask = nothing, accum = nothing, desc = nothing) -# eunion!(C, A, α, B, β, BinaryOp(op); mask, accum, desc) -# end -# -# function eunion(A, α, B, β, op::Function; mask = nothing, accum = nothing, desc = nothing) -# eunion(A, α, B, β, BinaryOp(op); mask, accum, desc) -# end - function Base.:+(A::GBArray, B::GBArray) eadd(A, B, +) end diff --git a/src/operations/extract.jl b/src/operations/extract.jl index 823a79a8..ae1f1ecd 100644 --- a/src/operations/extract.jl +++ b/src/operations/extract.jl @@ -54,7 +54,7 @@ Extract a submatrix or subvector from `A` into `C`. extract! function extract!( - C::GBMatrix, A::GBMatOrTranspose, I, J; + C::AbstractGBMatrix, A::GBMatOrTranspose, I, J; mask = nothing, accum = nothing, desc = nothing ) I, ni = idx(I) @@ -65,28 +65,28 @@ function extract!( desc = _handledescriptor(desc; in1 = A) I = decrement!(I) J = decrement!(J) - @wraperror LibGraphBLAS.GrB_Matrix_extract(C, mask, getaccum(accum, eltype(C)), parent(A), I, ni, J, nj, desc) + @wraperror LibGraphBLAS.GrB_Matrix_extract(gbpointer(C), mask, getaccum(accum, eltype(C)), gbpointer(parent(A)), I, ni, J, nj, desc) I isa Vector && increment!(I) J isa Vector && increment!(J) return C end function extract!( - C::GBMatrix, A::GBMatOrTranspose, ::Colon, J; + C::AbstractGBMatrix, A::GBMatOrTranspose, ::Colon, J; mask = nothing, accum = nothing, desc = nothing ) return extract!(C, A, ALL, J; mask, accum, desc) end function extract!( - C::GBMatrix, A::GBMatOrTranspose, I, ::Colon; + C::AbstractGBMatrix, A::GBMatOrTranspose, I, ::Colon; mask = nothing, accum = nothing, desc = nothing ) return extract!(C, A, I, ALL; mask, accum, desc) end function extract!( - C::GBMatrix, A::GBMatOrTranspose, ::Colon, ::Colon; + C::AbstractGBMatrix, A::GBMatOrTranspose, ::Colon, ::Colon; mask = nothing, accum = nothing, desc = nothing ) return extract!(C, A, ALL, ALL; mask, accum, desc) @@ -155,14 +155,14 @@ function Base.getindex( end function extract!( - w::GBVector, u::GBVector, I; + w::AbstractGBVector, u::AbstractGBVector, I; mask = nothing, accum = nothing, desc = nothing ) I, ni = idx(I) I = decrement!(I) desc = _handledescriptor(desc) mask === nothing && (mask = C_NULL) - @wraperror LibGraphBLAS.GrB_Matrix_extract(w, mask, getaccum(accum, eltype(w)), u, I, ni, UInt64[1], 1, desc) + @wraperror LibGraphBLAS.GrB_Matrix_extract(gbpointer(w), mask, getaccum(accum, eltype(w)), gbpointer(u), I, ni, UInt64[1], 1, desc) I isa Vector && increment!(I) return w end diff --git a/src/operations/kronecker.jl b/src/operations/kronecker.jl index 4b656ea8..b00b55c7 100644 --- a/src/operations/kronecker.jl +++ b/src/operations/kronecker.jl @@ -4,7 +4,7 @@ In-place version of [kron](@ref). """ function LinearAlgebra.kron!( - C::GBArray, + C::GBVecOrMat, A::GBArray, B::GBArray, op = *; @@ -16,15 +16,7 @@ function LinearAlgebra.kron!( desc = _handledescriptor(desc; in1=A, in2=B) op = BinaryOp(op)(eltype(A), eltype(B)) accum = getaccum(accum, eltype(C)) - if op isa TypedBinaryOperator - @wraperror LibGraphBLAS.GxB_kron(C, mask, accum, op, parent(A), parent(B), desc) - elseif op isa TypedMonoid - @wraperror LibGraphBLAS.GrB_Matrix_kronecker_Monoid(C, mask, accum, op, parent(A), parent(B), desc) - elseif op isa TypedSemiring - @wraperror LibGraphBLAS.GrB_Matrix_kronecker_Semiring(C, mask, accum, op, parent(A), parent(B), desc) - else - throw(ArgumentError("$op is not a valid monoid binary op or semiring.")) - end + @wraperror LibGraphBLAS.GxB_kron(gbpointer(C), mask, accum, op, gbpointer(parent(A)), gbpointer(parent(B)), desc) return C end @@ -55,7 +47,7 @@ function LinearAlgebra.kron( desc = nothing ) t = inferbinarytype(eltype(A), eltype(B), op) - C = GBMatrix{t}(size(A,1) * size(B, 1), size(A, 2) * size(B, 2)) + C = similar(A, t, (size(A, 1) * size(B, 1), size(A, 2) * size(B, 2)); fill = _promotefill(parent(A).fill, parent(B).fill)) kron!(C, A, B, op; mask, accum, desc) return C end \ No newline at end of file diff --git a/src/operations/map.jl b/src/operations/map.jl index 7d94d3ff..e829baee 100644 --- a/src/operations/map.jl +++ b/src/operations/map.jl @@ -6,7 +6,7 @@ function apply!( desc = _handledescriptor(desc; in1=A) op = UnaryOp(op)(eltype(A)) accum = getaccum(accum, eltype(C)) - @wraperror LibGraphBLAS.GrB_Matrix_apply(C, mask, accum, op, parent(A), desc) + @wraperror LibGraphBLAS.GrB_Matrix_apply(gbpointer(C), mask, accum, op, gbpointer(parent(A)), desc) return C end @@ -55,7 +55,7 @@ function apply!( desc = _handledescriptor(desc; in2=A) op = BinaryOp(op)(eltype(A), typeof(x)) accum = getaccum(accum, eltype(C)) - @wraperror LibGraphBLAS.GxB_Matrix_apply_BinaryOp1st(C, mask, accum, op, GBScalar(x), parent(A), desc) + @wraperror LibGraphBLAS.GxB_Matrix_apply_BinaryOp1st(gbpointer(C), mask, accum, op, GBScalar(x), gbpointer(parent(A)), desc) return C end @@ -82,7 +82,7 @@ function apply!( desc = _handledescriptor(desc; in1=A) op = BinaryOp(op)(eltype(A), typeof(x)) accum = getaccum(accum, eltype(C)) - @wraperror LibGraphBLAS.GxB_Matrix_apply_BinaryOp2nd(C, mask, accum, op, parent(A), GBScalar(x), desc) + @wraperror LibGraphBLAS.GxB_Matrix_apply_BinaryOp2nd(gbpointer(C), mask, accum, op, gbpointer(parent(A)), GBScalar(x), desc) return C end @@ -111,9 +111,9 @@ function Base.map!(f, A::GBArray{T}; mask = nothing, accum = nothing, desc = not apply!(f, C, A; mask, accum, desc) end -Base.:*(x, u::GBArray{T}; mask = nothing, accum = nothing, desc = nothing) where {T} = +Base.:*(x::V, u::GBArray{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; mask = nothing, accum = nothing, desc = nothing) where {T} = +Base.:*(u::GBArray{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) @@ -128,16 +128,13 @@ function mask!(C::GBArray, A::GBArray, mask::GBArray; structural = false, comple desc = Descriptor() structural && (desc.structural_mask=true) complement && (desc.complement_mask=true) + mask = mask isa Transpose || mask isa Adjoint ? copy(mask) : mask apply!(identity, C, A; mask, desc) return C end function mask!(A::GBArray, mask::GBArray; structural = false, complement = false) - desc = Descriptor() - structural && (desc.structural_mask=true) - complement && (desc.complement_mask=true) - apply!(identity, A, A; mask, desc) - return A + mask!(A, A, mask; structural, complement) end """ diff --git a/src/operations/mul.jl b/src/operations/mul.jl index 327396de..1f4e78d7 100644 --- a/src/operations/mul.jl +++ b/src/operations/mul.jl @@ -1,3 +1,7 @@ +# TODO: +# vxm shouldn't be done "magically" +# Instead it should be required that we have AdjTrans{<:GBVector} = AdjTrans{<:GBVector} * GBArray. + function LinearAlgebra.mul!( C::GBVecOrMat, A::GBArray, @@ -15,7 +19,7 @@ function LinearAlgebra.mul!( 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(C, mask, accum, op, parent(A), parent(B), desc) + @wraperror LibGraphBLAS.GrB_mxm(gbpointer(C), mask, accum, op, gbpointer(parent(A)), gbpointer(parent(B)), desc) return C end @@ -82,14 +86,15 @@ function mul( desc = nothing ) t = inferbinarytype(eltype(A), eltype(B), op) - if A isa GBMatOrTranspose && B isa GBVector - C = GBVector{t}(size(A, 1)) + fill = _promotefill(parent(A).fill, parent(B).fill) + if A isa GBMatOrTranspose && B isa AbstractGBVector + C = similar(A, t, size(A, 1); fill) elseif A isa GBVector && B isa GBMatOrTranspose - C = GBVector{t}(size(B, 2)) + C = similar(A, t, size(B, 2); fill) elseif A isa Transpose{<:Any, <:GBVector} && B isa GBVector - C = GBVector{t}(1) + C = similar(A, t, 1; fill) else - C = GBMatrix{t}(size(A, 1), size(B, 2)) + C = similar(A, t, (size(A, 1), size(B, 2)); fill) end mul!(C, A, B, op; mask, accum, desc) return C @@ -133,6 +138,28 @@ function Base.:*( return mul(A, B, (+, *); mask, accum, desc) end +# clear up some ambiguities: +function Base.:*( + A::AbstractGBVector{T}, + B::Transpose{T, <:AbstractGBVector{T}}; + mask = nothing, + accum = nothing, + desc = nothing +) where {T <: Real} + return mul(A, B, (+, *); mask, accum, desc) +end + +function Base.:*( + A::Transpose{T, <:AbstractGBVector{T}}, + B::AbstractGBVector{T}; + mask = nothing, + accum = nothing, + desc = nothing +) where {T <: Real} + return mul(A, B, (+, *); mask, accum, desc) +end + + function Base.:*( A::VecOrMat, B::GBArray; diff --git a/src/operations/operationutils.jl b/src/operations/operationutils.jl index 8bbced5a..d67e87ab 100644 --- a/src/operations/operationutils.jl +++ b/src/operations/operationutils.jl @@ -46,3 +46,16 @@ function ytype end Determine type of the output of a typed operator. """ function ztype end + +_promotefill(::Nothing, ::Nothing) = nothing +_promotefill(::Nothing, x) = nothing +_promotefill(x, ::Nothing) = nothing +_promotefill(::Missing, ::Missing) = missing +_promotefill(::Missing, x) = missing +_promotefill(x, ::Missing) = missing +# I'd prefer that this be nothing on x != y. But for type inference reasons this seems better. +# It's not a serious issue for several reasons. +# The first is that GrB methods don't know anything about fill, they don't care. +# The second is that it's free to setfill(A, nothing). Methods that are sensitive to this can enforce that. +# And third a future GBGraph type can manage this for the user. +_promotefill(x::X, y::Y) where {X, Y} = x == y ? (return promote_type(X, Y)(x)) : (return zero(promote_type(X, Y))) \ No newline at end of file diff --git a/src/operations/reduce.jl b/src/operations/reduce.jl index fc95aefb..6ee0bfbe 100644 --- a/src/operations/reduce.jl +++ b/src/operations/reduce.jl @@ -1,12 +1,14 @@ function reduce!( - op, w::GBVector, A::GBMatOrTranspose; + op, w::AbstractGBVector, A::GBMatOrTranspose; mask = nothing, accum = nothing, desc = nothing ) mask, accum = _handlenothings(mask, accum) desc = _handledescriptor(desc; in1=A) op = Monoid(op)(eltype(w)) accum = getaccum(accum, eltype(w)) - @wraperror LibGraphBLAS.GrB_Matrix_reduce_Monoid(Ptr{LibGraphBLAS.GrB_Vector}(w.p), mask, accum, op, parent(A), desc) + @wraperror LibGraphBLAS.GrB_Matrix_reduce_Monoid( + Ptr{LibGraphBLAS.GrB_Vector}(gbpointer(w)), mask, accum, op, gbpointer(parent(A)), desc + ) return w end @@ -26,12 +28,12 @@ function Base.reduce( end if dims == 2 && !(A isa GBVecOrTranspose) - w = GBVector{typeout}(size(A, 1)) + w = similar(A, typeout, size(A, 1)) reduce!(op, w, A; desc, accum, mask) return w elseif dims == 1 && !(A isa GBVecOrTranspose) desc.transpose_input1 = true - w = GBVector{typeout}(size(A, 2)) + 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 @@ -45,7 +47,7 @@ function Base.reduce( op = Monoid(op)(typec) desc = _handledescriptor(desc; in1=A) accum = getaccum(accum, typec) - @wraperror LibGraphBLAS.GrB_Matrix_reduce_Monoid_Scalar(c, accum, op, parent(A), desc) + @wraperror LibGraphBLAS.GrB_Matrix_reduce_Monoid_Scalar(c, accum, op, gbpointer(parent(A)), desc) return c[] end end diff --git a/src/operations/resize.jl b/src/operations/resize.jl index 06caf203..0538a519 100644 --- a/src/operations/resize.jl +++ b/src/operations/resize.jl @@ -1,9 +1,9 @@ -function Base.resize!(A::GBMatrix, nrows::Integer, ncols::Integer) - @wraperror LibGraphBLAS.GrB_Matrix_resize(A, LibGraphBLAS.GrB_Index(nrows), LibGraphBLAS.GrB_Index(ncols)) +function Base.resize!(A::AbstractGBMatrix, nrows::Integer, ncols::Integer) + @wraperror LibGraphBLAS.GrB_Matrix_resize(gbpointer(A), LibGraphBLAS.GrB_Index(nrows), LibGraphBLAS.GrB_Index(ncols)) return A end -function Base.resize!(v::GBVector, nrows::Integer) - @wraperror LibGraphBLAS.GrB_Matrix_resize(v, LibGraphBLAS.GrB_Index(nrows), 1) +function Base.resize!(v::AbstractGBVector, nrows::Integer) + @wraperror LibGraphBLAS.GrB_Matrix_resize(gbpointer(v), LibGraphBLAS.GrB_Index(nrows), 1) return v end diff --git a/src/operations/select.jl b/src/operations/select.jl index b076d0ce..5a4b0c51 100644 --- a/src/operations/select.jl +++ b/src/operations/select.jl @@ -1,6 +1,8 @@ +# TODO: update to modern op system. + "In place version of `select`." function select!( - op::SelectUnion, + op, C::GBVecOrMat, A::GBArray, thunk = nothing; @@ -8,6 +10,7 @@ function select!( accum = nothing, desc = nothing ) + op = SelectOp(op) mask, accum = _handlenothings(mask, accum) desc = _handledescriptor(desc; in1=A) thunk === nothing && (thunk = C_NULL) @@ -15,22 +18,10 @@ function select!( if thunk isa Number thunk = GBScalar(thunk) end - @wraperror LibGraphBLAS.GxB_Matrix_select(C, mask, accum, op, parent(A), thunk, desc) + @wraperror LibGraphBLAS.GxB_Matrix_select(C, mask, accum, op, gbpointer(parent(A)), thunk, desc) return C end -function select!( - op, - C::GBVecOrMat, - A::GBArray, - thunk = nothing; - mask = nothing, - accum = nothing, - desc = nothing -) - return select!(SelectOp(op), C, A, thunk; mask, accum, desc) -end - function select!(op, A::GBArray, thunk = nothing; mask = nothing, accum = nothing, desc = nothing) return select!(op, A, A, thunk; mask, accum, desc) end @@ -60,35 +51,19 @@ Some SelectOps or functions may require an additional argument `thunk`, for use - `GBArray`: The output matrix whose `eltype` is determined by `A` and `op`. """ function select( - op::SelectUnion, + op, A::GBArray, - thunk::Union{GBScalar, Nothing, valid_union} = nothing; + thunk = nothing; mask = nothing, accum = nothing, desc = nothing ) + op = SelectOp(op) mask, accum = _handlenothings(mask, accum) C = similar(A) select!(op, C, A, thunk; accum, mask, desc) return C end -function select( - op::Function, A::GBArray, thunk; - mask = nothing, accum = nothing, desc = nothing -) - select(SelectOp(op), A, thunk; mask, accum, desc) -end - -function select( - op::Function, - A::GBArray, - thunk::Union{GBScalar, Nothing, Number} = nothing; - mask = nothing, - accum = nothing, - desc = nothing -) - return select(SelectOp(op), A, thunk; mask, accum, desc) -end LinearAlgebra.tril(A::GBArray, k::Integer = 0) = select(tril, A, k) LinearAlgebra.triu(A::GBArray, k::Integer = 0) = select(triu, A, k) diff --git a/src/operations/sort.jl b/src/operations/sort.jl index b169593b..febd0474 100644 --- a/src/operations/sort.jl +++ b/src/operations/sort.jl @@ -24,7 +24,7 @@ function Base.sort!( end desc = _handledescriptor(desc; in1=A) desc.transpose_input1 = transpose - @wraperror LibGraphBLAS.GxB_Matrix_sort(C, P, op, parent(A), desc) + @wraperror LibGraphBLAS.GxB_Matrix_sort(C, P, op, gbpointer(parent(A)), desc) return C end diff --git a/src/operations/transpose.jl b/src/operations/transpose.jl index d4704522..54112bb7 100644 --- a/src/operations/transpose.jl +++ b/src/operations/transpose.jl @@ -15,16 +15,15 @@ Eagerly evaluated matrix transpose, storing the output in `C`. - `desc::Union{Nothing, Descriptor} = DEFAULTDESC` """ function gbtranspose!( - C::GBVecOrMat, A::GBArray; + C::AbstractGBArray, A::GBArray; mask = nothing, accum = nothing, desc = nothing ) mask, accum = _handlenothings(mask, accum) desc = _handledescriptor(desc; in1=A) accum = getaccum(accum, eltype(C)) - @wraperror LibGraphBLAS.GrB_transpose(C, mask, accum, parent(A), desc) + @wraperror LibGraphBLAS.GrB_transpose(gbpointer(C), mask, accum, gbpointer(parent(A)), desc) return C end - """ gbtranspose(A::GBMatrix; kwargs...)::GBMatrix @@ -39,10 +38,7 @@ Eagerly evaluated matrix transpose which returns the transposed matrix. # Returns - `C::GBMatrix`: output matrix. """ -function gbtranspose( - A::GBArray; - mask = nothing, accum = nothing, desc = nothing -) +function gbtranspose(A::GBArray; mask = nothing, accum = nothing, desc = nothing) C = similar(A, size(A,2), size(A, 1)) gbtranspose!(C, A; mask, accum, desc) return C @@ -53,25 +49,21 @@ function LinearAlgebra.transpose(A::GBArray) end function Base.copy!( - C::GBMatrix, A::LinearAlgebra.Transpose{<:Any, <:GBArray}; + C::GBVecOrMat, A::LinearAlgebra.Transpose{<:Any, <:GBVecOrMat}; mask = nothing, accum = nothing, desc = nothing ) return gbtranspose!(C, A.parent; mask, accum, desc) end - - function Base.copy( - A::LinearAlgebra.Transpose{<:Any, <:GBArray}; + A::LinearAlgebra.Transpose{<:Any, <:GBVecOrMat}; mask = nothing, accum = nothing, desc = nothing ) return gbtranspose(parent(A); mask, accum, desc) end #This is ok per the GraphBLAS Slack channel. Should change its effect on Complex input. -LinearAlgebra.adjoint(A::GBMatrix) = transpose(A) - -LinearAlgebra.adjoint(v::GBVector) = transpose(v) +LinearAlgebra.adjoint(A::GBVecOrMat) = transpose(A) #arrrrgh, type piracy. LinearAlgebra.transpose(::Nothing) = nothing diff --git a/src/operators/selectops.jl b/src/operators/selectops.jl index 13dac5d6..f90c61a3 100644 --- a/src/operators/selectops.jl +++ b/src/operators/selectops.jl @@ -11,6 +11,8 @@ mutable struct SelectOp <: AbstractSelectOp end end +SelectOp(op::SelectOp) = op + const SelectUnion = Union{AbstractSelectOp, LibGraphBLAS.GxB_SelectOp} Base.unsafe_convert(::Type{LibGraphBLAS.GxB_SelectOp}, selectop::SelectOp) = selectop.p diff --git a/src/options.jl b/src/options.jl index 2a49b691..ce89bff9 100644 --- a/src/options.jl +++ b/src/options.jl @@ -57,7 +57,7 @@ function GxB_Global_Option_set(field, value) end end -function GxB_Matrix_Option_get(A, field) +function GxB_Matrix_Option_get(A::AbstractGBArray, field) if field ∈ [GxB_HYPER_SWITCH, GxB_BITMAP_SWITCH] T = Cdouble elseif field ∈ [GxB_FORMAT, GxB_SPARSITY_STATUS, GxB_SPARSITY_CONTROL] @@ -68,20 +68,20 @@ function GxB_Matrix_Option_get(A, field) (:GxB_Matrix_Option_get, libgraphblas), Cvoid, (LibGraphBLAS.GrB_Matrix, UInt32, Ptr{Cvoid}), - A, + gbpointer(A), field, v ) return v[] end -function GxB_Matrix_Option_set(A, field, value) +function GxB_Matrix_Option_set(A::AbstractGBArray, field, value) if field ∈ [GxB_HYPER_SWITCH, GxB_BITMAP_SWITCH] ccall( (:GxB_Matrix_Option_set, libgraphblas), Cvoid, (LibGraphBLAS.GrB_Matrix, UInt32, Cdouble), - A, + gbpointer(A), field, value ) @@ -90,58 +90,7 @@ function GxB_Matrix_Option_set(A, field, value) (:GxB_Matrix_Option_set, libgraphblas), Cvoid, (LibGraphBLAS.GrB_Matrix, UInt32, UInt32), - A, - field, - value - ) - end -end - -function GxB_Vector_Option_get(A, field) - if field ∈ [GxB_HYPER_SWITCH, GxB_BITMAP_SWITCH] - T = Cdouble - elseif field ∈ [GxB_FORMAT] - T = UInt32 - elseif field ∈ [GxB_SPARSITY_STATUS, GxB_SPARSITY_CONTROL] - T = Cint - end - v = Ref{T}() - ccall( - (:GxB_Vector_Option_get, libgraphblas), - Cvoid, - (LibGraphBLAS.GrB_Vector, UInt32, Ptr{Cvoid}), - A, - field, - v - ) - return v[] -end - -function GxB_Vector_Option_set(A, field, value) - if field ∈ [GxB_HYPER_SWITCH, GxB_BITMAP_SWITCH] - ccall( - (:GxB_Vector_Option_set, libgraphblas), - Cvoid, - (LibGraphBLAS.GrB_Vector, UInt32, Cdouble), - A, - field, - value - ) - elseif field ∈ [GxB_FORMAT] - ccall( - (:GxB_Vector_Option_set, libgraphblas), - Cvoid, - (LibGraphBLAS.GrB_Vector, UInt32, UInt32), - A, - field, - value - ) - elseif field ∈ [GxB_SPARSITY_CONTROL] - ccall( - (:GxB_Vector_Option_set, libgraphblas), - Cvoid, - (LibGraphBLAS.GrB_Vector, UInt32, Cint), - A, + gbpointer(A), field, value ) @@ -160,30 +109,18 @@ function gbget(option) return GxB_Global_Option_get(option) end -function gbset(A::GBMatrix, option, value) +function gbset(A::AbstractGBArray, option, value) option = option_toconst(option) value = option_toconst(value) GxB_Matrix_Option_set(A, option, value) return nothing end -function gbget(A::GBMatrix, option) +function gbget(A::AbstractGBArray, option) option = option_toconst(option) return GxB_Matrix_Option_get(A, option) end -function gbset(A::GBVector, option, value) - option = option_toconst(option) - value = option_toconst(value) - GxB_Matrix_Option_set(A.p, option, value) - return nothing -end - -function gbget(A::GBVector, option) - option = option_toconst(option) - return GxB_Matrix_Option_get(A.p, option) -end - function format(A::GBVecOrMat) t = gbget(A, SPARSITY_STATUS) f = gbget(A, FORMAT) diff --git a/src/pack.jl b/src/pack.jl index 0d566e24..bd1d96bd 100644 --- a/src/pack.jl +++ b/src/pack.jl @@ -1,10 +1,10 @@ -function _packdensematrix!(A::GBVecOrMat{T}, M::DenseVecOrMat; desc = nothing) where {T} +function _packdensematrix!(A::AbstractGBArray{T}, M::DenseVecOrMat; desc = nothing) where {T} desc = _handledescriptor(desc) Csize = length(A) * sizeof(T) values = Ref{Ptr{Cvoid}}(pointer(M)) isuniform = false @wraperror LibGraphBLAS.GxB_Matrix_pack_FullC( - A.p, + gbpointer(A), values, Csize, isuniform, @@ -23,15 +23,15 @@ function _packcscmatrix!( rowidxsize = length(rowidx) * sizeof(LibGraphBLAS.GrB_Index), valsize = length(values) * sizeof(T) ) where {T, Ti} - colptr .-= 1 - rowidx .-= 1 + decrement!(colptr) + decrement!(rowidx) colptr = Ref{Ptr{LibGraphBLAS.GrB_Index}}(pointer(colptr)) rowidx = Ref{Ptr{LibGraphBLAS.GrB_Index}}(pointer(rowidx)) values = Ref{Ptr{Cvoid}}(pointer(values)) desc = _handledescriptor(desc) @wraperror LibGraphBLAS.GxB_Matrix_pack_CSC( - A, + gbpointer(A), colptr, rowidx, values, @@ -55,15 +55,15 @@ function _packcsrmatrix!( colidxsize = length(colidx) * sizeof(LibGraphBLAS.GrB_Index), valsize = length(values) * sizeof(T) ) where {T, Ti} - rowptr .-= 1 - colidx .-= 1 + decrement!(rowptr) + decrement!(colidx) rowptr = Ref{Ptr{LibGraphBLAS.GrB_Index}}(pointer(rowptr)) colidx = Ref{Ptr{LibGraphBLAS.GrB_Index}}(pointer(colidx)) values = Ref{Ptr{Cvoid}}(pointer(values)) desc = _handledescriptor(desc) @wraperror LibGraphBLAS.GxB_Matrix_pack_CSC( - A, + gbpointer(A), rowptr, colidx, values, @@ -77,6 +77,6 @@ function _packcsrmatrix!( return A end -function _makeshallow!(A::GBVecOrMat) - ccall((:GB_make_shallow, libgraphblas), Cvoid, (LibGraphBLAS.GrB_Matrix,), A) +function _makeshallow!(A::AbstractGBArray) + ccall((:GB_make_shallow, libgraphblas), Cvoid, (LibGraphBLAS.GrB_Matrix,), gbpointer(A)) end \ No newline at end of file diff --git a/src/print.jl b/src/print.jl index 63d648bd..cc489866 100644 --- a/src/print.jl +++ b/src/print.jl @@ -19,17 +19,17 @@ function gxbstring(x, name = "", level::LibGraphBLAS.GxB_Print_Level = LibGraphB @wraperror LibGraphBLAS.GxB_Semiring_fprint(x, name, level, cf) elseif x isa AbstractDescriptor @wraperror LibGraphBLAS.GxB_Descriptor_fprint(x, name, level, cf) - elseif x isa GBVector + elseif x isa AbstractGBVector if level == LibGraphBLAS.GxB_SUMMARY - @wraperror LibGraphBLAS.GxB_Matrix_fprint(x, name, LibGraphBLAS.GxB_SHORT, cf) + @wraperror LibGraphBLAS.GxB_Matrix_fprint(gbpointer(x), name, LibGraphBLAS.GxB_SHORT, cf) else - @wraperror LibGraphBLAS.GxB_Matrix_fprint(x, name, level, cf) + @wraperror LibGraphBLAS.GxB_Matrix_fprint(gbpointer(x), name, level, cf) end - elseif x isa GBMatrix + elseif x isa AbstractGBMatrix if level == LibGraphBLAS.GxB_SUMMARY - @wraperror LibGraphBLAS.GxB_Matrix_fprint(x, name, LibGraphBLAS.GxB_SHORT, cf) + @wraperror LibGraphBLAS.GxB_Matrix_fprint(gbpointer(x), name, LibGraphBLAS.GxB_SHORT, cf) else - @wraperror LibGraphBLAS.GxB_Vector_fprint(x, name, level, cf) + @wraperror LibGraphBLAS.GxB_Vector_fprint(gbpointer(x), name, level, cf) end elseif x isa GBScalar @wraperror LibGraphBLAS.GxB_Scalar_fprint(x, name, LibGraphBLAS.GxB_SHORT, cf) diff --git a/src/random.jl b/src/random.jl index 06cc39aa..81865338 100644 --- a/src/random.jl +++ b/src/random.jl @@ -44,7 +44,7 @@ function gbrand( if !(type <: Complex) hermitian = false end - + # TODO: switch from A[i, j] = x, to COO->build for _ ∈ 1:round(Int64, nrows * ncols * density) i = rand(rng, 1:nrows) j = rand(rng, 1:ncols) diff --git a/src/types.jl b/src/types.jl index a086975b..b5ffecb8 100644 --- a/src/types.jl +++ b/src/types.jl @@ -257,19 +257,9 @@ compressed sparse vector. See also: [`GBMatrix`](@ref). """ -mutable struct GBVector{T} <: AbstractSparseArray{T, UInt64, 1} - p::LibGraphBLAS.GrB_Matrix - function GBVector{T}(p::LibGraphBLAS.GrB_Matrix; aliased=false) where {T} - v = new(p) - function f(vector) - @wraperror LibGraphBLAS.GrB_Matrix_free(Ref(vector.p)) - end - if aliased - return v - else - return finalizer(f, v) - end - end +mutable struct GBVector{T, F} <: AbstractGBVector{T, F} + p::Ref{LibGraphBLAS.GrB_Matrix} # a GBVector is a GBMatrix internally. + fill::F end """ @@ -286,19 +276,12 @@ row or column orientation: The storage type is automatically determined by the library. """ -mutable struct GBMatrix{T} <: AbstractSparseArray{T, UInt64, 2} - p::LibGraphBLAS.GrB_Matrix - # NOTE WELL: The alias kwarg IS NOT for public use. - # It is used in very few places to convert a GBVector to a GBMatrix internally. - function GBMatrix{T}(p::LibGraphBLAS.GrB_Matrix; aliased=false) where {T} - A = new(p) - function f(matrix) - @wraperror LibGraphBLAS.GrB_Matrix_free(Ref(matrix.p)) - end - if aliased - return A - else - return finalizer(f, A) - end - end +mutable struct GBMatrix{T, F} <: AbstractGBMatrix{T, F} + p::Ref{LibGraphBLAS.GrB_Matrix} + fill::F end + +# Most likely this will be the case for all AbstractGBArray. +# However, if one (for some reason) wraps another AbstractGBArray +# this should be overloaded. +gbpointer(A::AbstractGBArray) = A.p[] \ No newline at end of file diff --git a/src/unpack.jl b/src/unpack.jl index a4b4790d..4b3caa84 100644 --- a/src/unpack.jl +++ b/src/unpack.jl @@ -1,10 +1,10 @@ -function _unpackdensematrix!(A::GBVecOrMat{T}; desc = nothing) where {T} +function _unpackdensematrix!(A::AbstractGBArray{T}; desc = nothing) where {T} desc = _handledescriptor(desc) Csize = Ref{LibGraphBLAS.GrB_Index}(length(A) * sizeof(T)) values = Ref{Ptr{Cvoid}}(Ptr{T}()) isiso = Ref{Bool}(false) @wraperror LibGraphBLAS.GxB_Matrix_unpack_FullC( - A.p, + gbpointer(A), values, Csize, isiso, @@ -13,7 +13,7 @@ function _unpackdensematrix!(A::GBVecOrMat{T}; desc = nothing) where {T} return unsafe_wrap(Array{T}, Ptr{T}(values[]), size(A)) end -function _unpackcscmatrix!(A::GBVecOrMat{T}; desc = nothing) where {T} +function _unpackcscmatrix!(A::AbstractGBArray{T}; desc = nothing) where {T} desc = _handledescriptor(desc) colptr = Ref{Ptr{LibGraphBLAS.GrB_Index}}() rowidx = Ref{Ptr{LibGraphBLAS.GrB_Index}}() @@ -25,7 +25,7 @@ function _unpackcscmatrix!(A::GBVecOrMat{T}; desc = nothing) where {T} isjumbled = C_NULL nnonzeros = nnz(A) @wraperror LibGraphBLAS.GxB_Matrix_unpack_CSC( - A.p, + gbpointer(A), colptr, rowidx, values, @@ -49,8 +49,8 @@ function _unpackcscmatrix!(A::GBVecOrMat{T}; desc = nothing) where {T} vals = unsafe_wrap(Array{T}, Ptr{T}(values[]), nnonzeros) valsize = valsize[] end - colptr .+= 1 - rowidx .+= 1 + increment!(colptr) + increment!(rowidx) return colptr, colptrsize[], rowidx, @@ -71,7 +71,7 @@ function _unpackcsrmatrix!(A::GBVecOrMat{T}; desc = nothing) where {T} isjumbled = C_NULL nnonzeros = nnz(A) @wraperror LibGraphBLAS.GxB_Matrix_unpack_CSR( - A.p, + gbpointer(A), rowptr, colidx, values, @@ -98,9 +98,9 @@ function _unpackcsrmatrix!(A::GBVecOrMat{T}; desc = nothing) where {T} vals = unsafe_wrap(Array{T}, Ptr{T}(values[]), nnonzeros) valsize = valsize[] end + increment!(rowptr) + increment!(colidx) - rowptr .+= 1 - colidx .+= 1 return rowptr, rowptrsize[], colidx, diff --git a/src/vector.jl b/src/vector.jl index aabc4d7a..63ae575c 100644 --- a/src/vector.jl +++ b/src/vector.jl @@ -1,30 +1,34 @@ # Constructors: ############### """ - GBVector{T}(n = LibGraphBLAS.GxB_INDEX_MAX) + GBVector{T}(n; fill = nothing) """ -function GBVector{T}(n = LibGraphBLAS.GxB_INDEX_MAX) where {T} +function GBVector{T}(n; fill::F = nothing) where {T, F} m = Ref{LibGraphBLAS.GrB_Matrix}() @wraperror LibGraphBLAS.GrB_Matrix_new(m, gbtype(T), n, 1) - v = GBVector{T}(m[]) + v = GBVector{T, F}(finalizer(m) do ref + @wraperror LibGraphBLAS.GrB_Matrix_free(ref) + end, fill) gbset(v, FORMAT, BYCOL) return v end -GBVector{T}(dims::Dims{1}) where {T} = GBVector{T}(dims...) -GBVector{T}(nrows::Base.OneTo) where {T} = - GBVector{T}(nrows.stop) -GBVector{T}(nrows::Tuple{Base.OneTo,}) where {T} = GBVector{T}(first(nrows)) +GBVector{T}(dims::Dims{1}; fill = nothing) where {T} = GBVector{T}(dims...; fill) +GBVector{T}(nrows::Base.OneTo; fill = nothing) where {T} = + GBVector{T}(nrows.stop; fill) +GBVector{T}(nrows::Tuple{Base.OneTo,}; fill = nothing) where {T} = GBVector{T}(first(nrows); fill) """ GBVector(I::AbstractVector, X::AbstractVector{T}) Create a GBVector from a vector of indices `I` and a vector of values `X`. """ -function GBVector(I::AbstractVector{U}, X::AbstractVector{T}; combine = +, nrows = maximum(I)) where {U<:Integer, T} - x = GBVector{T}(nrows) - build(x, I, X, combine = combine) - return x +function GBVector(I::AbstractVector{U}, X::AbstractVector{T}; combine = +, nrows = maximum(I), fill = nothing) where {U<:Integer, T} + I isa Vector || (I = collect(I)) + X isa Vector || (X = collect(X)) + v = GBVector{T}(nrows; fill) + build(v, I, X; combine) + return v end #iso valued constructors. @@ -36,8 +40,8 @@ The resulting vector is "iso-valued" such that it only stores `x` once rather th each index. """ function GBVector(I::AbstractVector{U}, x::T; - nrows = maximum(I)) where {U<:Integer, T} - A = GBVector{T}(nrows) + nrows = maximum(I), fill = nothing) where {U<:Integer, T} + A = GBVector{T}(nrows; fill) build(A, I, x) return A end @@ -49,73 +53,68 @@ Create an `n` length dense GBVector `v` such that M[I[k]] = x. The resulting vector is "iso-valued" such that it only stores `x` once rather than once for each index. """ -function GBVector(n::Integer, x::T) where {T} - v = GBVector{T}(n) +function GBVector(n::Integer, x::T; fill = nothing) where {T} + v = GBVector{T}(n; fill) v[:] = x return v end # Some Base and basic SparseArrays/LinearAlgebra functions: ########################################################### -Base.unsafe_convert(::Type{LibGraphBLAS.GrB_Matrix}, v::GBVector) = v.p +Base.unsafe_convert(::Type{LibGraphBLAS.GrB_Matrix}, v::GBVector) = v.p[] -function Base.copy(A::GBVector{T}) where {T} +function Base.copy(A::GBVector{T, F}) where {T, F} C = Ref{LibGraphBLAS.GrB_Matrix}() - LibGraphBLAS.GrB_Matrix_dup(C, A) - return GBVector{T}(C[]) -end - -function Base.size(v::GBVector) - nrows = Ref{LibGraphBLAS.GrB_Index}() - @wraperror LibGraphBLAS.GrB_Matrix_nrows(nrows, v) - return (Int64(nrows[]),) + LibGraphBLAS.GrB_Matrix_dup(C, gbpointer(A)) + return GBVector{T, F}(C[], A.fill) end -Base.eltype(::Type{GBVector{T}}) where{T} = T +# because of the fill kwarg we have to redo a lot of the Base.similar dispatch stack. function Base.similar( - ::GBVector{T}, ::Type{TNew}, - dims::Dims{1} -) where {T, TNew} - return GBVector{TNew}(dims...) + v::GBVectorOrTranspose{T}, ::Type{TNew} = T, + dims::Tuple{Int64, Vararg{Int64, N}} = size(v); fill = parent(v).fill +) where {T, TNew, N} + if dims isa Dims{1} + return GBVector{TNew}(dims...; fill) + else + return GBMatrix{TNew}(dims...; fill) + end +end + +function Base.similar(v::GBVectorOrTranspose{T}, dims::Tuple; fill = v.fill) where {T} + return similar(v, T, dims; fill) end function Base.similar( - ::GBVector{T}, ::Type{TNew}, - dims::Dims{2} + v::GBVectorOrTranspose{T}, ::Type{TNew}, + dims::Integer; fill = parent(v).fill ) where {T, TNew} - return GBMatrix{TNew}(dims...) + return similar(v, TNew, (dims,); fill) end -function Base.deleteat!(v::GBVector, i) - @wraperror LibGraphBLAS.GrB_Matrix_removeElement(v, decrement!(i), 1) - return v +function Base.similar( + v::GBVectorOrTranspose{T}, + dims::Integer; fill = parent(v).fill +) where {T} + return similar(v, (dims,); fill) end -function Base.resize!(v::GBVector, n) - @wraperror LibGraphBLAS.GrB_Matrix_resize(v, n, 1) - return v +function Base.similar( + v::GBVectorOrTranspose{T}, ::Type{TNew}, + dim1::Integer, dim2::Integer; fill = parent(v).fill +) where {T, TNew} + return similar(v, TNew, (dim1, dim2); fill) end -function LinearAlgebra.diag(A::GBMatOrTranspose{T}, k::Integer = 0; desc = nothing) where {T} - m, n = size(A) - if !(k in -m:n) - s = 0 - elseif k >= 0 - s = min(m, n - k) - else - s = min(m + k, n) - end - v = GBVector{T}(s) - desc = _handledescriptor(desc; in1=A) - if A isa Transpose - k = -k - end - @wraperror LibGraphBLAS.GxB_Vector_diag(LibGraphBLAS.GrB_Vector(v.p), parent(A), k, desc) - return v +function Base.similar( + v::GBVectorOrTranspose{T}, dim1::Integer, dim2::Integer; fill = parent(v).fill +) where {T} + return similar(v, (dim1, dim2); fill) end #We need these until I can get a SparseArrays.nonzeros implementation +# TODO: REMOVE function Base.show(io::IO, ::MIME"text/plain", v::GBVector) gxbprint(io, v) end @@ -128,100 +127,6 @@ function Base.show(io::IOContext, v::GBVector) gxbprint(io, v) end -# Type dependent functions build, setindex, getindex, and findnz: -for T ∈ valid_vec - if T ∈ gxb_vec - prefix = :GxB - else - prefix = :GrB - end - - # Build functions - func = Symbol(prefix, :_Matrix_build_, suffix(T)) - @eval begin - function build(v::GBVector{$T}, I::Vector, X::Vector{$T}; combine = +) - nnz(v) == 0 || throw(OutputNotEmptyError("Cannot build vector with existing elements")) - length(X) == length(I) || DimensionMismatch("I and X must have the same length") - combine = BinaryOp(combine)($T) - decrement!(I) - @wraperror LibGraphBLAS.$func( - Ptr{LibGraphBLAS.GrB_Vector}(v.p), - Vector{LibGraphBLAS.GrB_Index}(I), - zeros(LibGraphBLAS.GrB_Index, length(I)), - X, - length(X), - combine - ) - increment!(I) - end - end - # Setindex functions - func = Symbol(prefix, :_Matrix_setElement_, suffix(T)) - @eval begin - function Base.setindex!(v::GBVector{$T}, x, i::Integer) - x = convert($T, x) - return LibGraphBLAS.$func(v, x, LibGraphBLAS.GrB_Index(decrement!(i)), 0) - end - end - # Getindex functions - func = Symbol(prefix, :_Matrix_extractElement_, suffix(T)) - @eval begin - function Base.getindex(v::GBVector{$T}, i::Integer) - x = Ref{$T}() - result = LibGraphBLAS.$func(x, v, LibGraphBLAS.GrB_Index(decrement!(i)), 0) - if result == LibGraphBLAS.GrB_SUCCESS - return x[] - elseif result == LibGraphBLAS.GrB_NO_VALUE - return nothing - else - @wraperror result - end - end - end - # findnz functions - func = Symbol(prefix, :_Matrix_extractTuples_, suffix(T)) - @eval begin - function SparseArrays.findnz(v::GBVector{$T}) - nvals = Ref{LibGraphBLAS.GrB_Index}(nnz(v)) - I = Vector{LibGraphBLAS.GrB_Index}(undef, nvals[]) - X = Vector{$T}(undef, nvals[]) - @wraperror LibGraphBLAS.$func(I, C_NULL, X, nvals, v) - nvals[] == length(I) == length(X) || throw(DimensionMismatch("length(I) != length(X)")) - return increment!(I), X - end - function SparseArrays.nonzeros(v::GBVector{$T}) - nvals = Ref{LibGraphBLAS.GrB_Index}(nnz(v)) - X = Vector{$T}(undef, nvals[]) - @wraperror LibGraphBLAS.$func(C_NULL, C_NULL, X, nvals, v) - nvals[] == length(X) || throw(DimensionMismatch("")) - return X - end - function SparseArrays.nonzeroinds(v::GBVector{$T}) - nvals = Ref{LibGraphBLAS.GrB_Index}(nnz(v)) - I = Vector{LibGraphBLAS.GrB_Index}(undef, nvals[]) - wait(v) - @wraperror LibGraphBLAS.$func(I, C_NULL, C_NULL, nvals, v) - nvals[] == length(I) || throw(DimensionMismatch("")) - return increment!(I) - end - end -end - -function build(v::GBVector{T}, I::Vector, x::T) where {T} - nnz(v) == 0 || throw(OutputNotEmptyError("Cannot build vector with existing elements")) - x = GBScalar(x) - decrement!(I) - @wraperror LibGraphBLAS.GxB_Matrix_build_Scalar( - v, - Vector{LibGraphBLAS.GrB_Index}(I), - zeros(LibGraphBLAS.GrB_Index, length(I)), - x, - length(I) - ) - increment!(I) - return v -end - # Indexing functions: ##################### @@ -242,43 +147,3 @@ function Base.getindex( ) return extract(u, i; mask, accum, desc) end - -""" - subassign(w::GBVector, u::GBVector, I; kwargs...)::GBVector - -Assign a subvector of `w` to `u`. Return `u`. Equivalent to the matrix definition. -""" -function subassign!(w::GBVector{T}, u, I; mask = nothing, accum = nothing, desc = nothing) where {T} - return subassign!(GBMatrix{T}(w.p; aliased=true), u, I, UInt64[1]; mask, accum, desc) -end - -""" - assign(w::GBVector, u::GBVector, I; kwargs...)::GBVector - -Assign a subvector of `w` to `u`. Return `u`. Equivalent to the matrix definition. -""" -function assign!(w::GBVector{T}, u, I; mask = nothing, accum = nothing, desc = nothing) where {T} - return assign!(GBMatrix{T}(w.p; aliased=true), u, I, UInt64[1]; mask, accum, desc) -end - -function Base.setindex!( - u::GBVector, x, ::Colon; - mask = nothing, accum = nothing, desc = nothing -) - subassign!(u, x, ALL; mask, accum, desc) - return nothing -end -# silly overload to help a bit with broadcasting. -function Base.setindex!( - u::GBVector, x, I::Union{Vector, UnitRange, StepRange, Colon}, ::Colon; - mask = nothing, accum = nothing, desc = nothing -) - Base.setindex!(u, x, I; mask, accum, desc) -end -function Base.setindex!( - u::GBVector, x, I::Union{Vector, UnitRange, StepRange}; - mask = nothing, accum = nothing, desc = nothing -) - subassign!(u, x, I; mask, accum, desc) - return nothing -end diff --git a/src/wait.jl b/src/wait.jl index 763c2dd3..c20167b1 100644 --- a/src/wait.jl +++ b/src/wait.jl @@ -1,6 +1,6 @@ -function Base.wait(A::GBArray) +function Base.wait(A::AbstractGBArray) waitmode = LibGraphBLAS.GrB_MATERIALIZE - @wraperror LibGraphBLAS.GrB_Matrix_wait(A, waitmode) + @wraperror LibGraphBLAS.GrB_Matrix_wait(gbpointer(A), waitmode) return nothing end diff --git a/test/operatorutils.jl b/test/operatorutils.jl index e40ffb0e..229faf3e 100644 --- a/test/operatorutils.jl +++ b/test/operatorutils.jl @@ -1,8 +1,8 @@ @testset "Operator Utilities" begin @test SuiteSparseGraphBLAS.optype(Float64, Float32) == Float64 @test SuiteSparseGraphBLAS.optype(UInt32, UInt64) == UInt64 - A = GBMatrix{Float64}() - B = GBMatrix{Int32}() + A = GBMatrix{Float64}(3, 3) + B = GBMatrix{Int32}(3, 3) @test SuiteSparseGraphBLAS.optype(A, B) == Float64 @test SuiteSparseGraphBLAS.symtotype(:nB) == SuiteSparseGraphBLAS.nBtypes diff --git a/test/runtests.jl b/test/runtests.jl index 97b928ee..1949261a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -90,7 +90,7 @@ println("Testing SuiteSparseGraphBLAS.jl") include_test("chainrules/selectrules.jl") include_test("chainrules/constructorrules.jl") include_test("chainrules/maprules.jl") - include_test("spmgb/sparsemat.jl") - include_test("spmgb/higherorderfns.jl") + # include_test("spmgb/sparsemat.jl") + # include_test("spmgb/higherorderfns.jl") end