Skip to content

Commit 8274567

Browse files
rayegunrayegun
committed
reduce cholesky manual labor, still a little buggy
1 parent 84b21d7 commit 8274567

File tree

15 files changed

+1004
-141
lines changed

15 files changed

+1004
-141
lines changed

‎src/abstractgbarray.jl‎

Lines changed: 31 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -28,7 +28,38 @@ function SparseArrays.nnz(A::GBArrayOrTranspose)
2828
return Int64(nvals[])
2929
end
3030

31+
function strip_parameters end
32+
promote_storage(::S, ::S) where {S <: StorageOrders.StorageOrder} = S()
33+
promote_storage(::S1, ::S2) where {S1 <: StorageOrders.StorageOrder, S2 <:StorageOrders.StorageOrder} =
34+
StorageOrders.RuntimeOrder()
35+
3136
# Base functions:
37+
# default to GBMatrix:
38+
Base.promote_rule(::Type{<:AbstractGBMatrix{T, F}}, ::Type{<:AbstractGBMatrix{T2, F2}}) where {T, F, T2, F2} =
39+
GBMatrix{promote_type(T, T2), promote_type(F, F2)}
40+
Base.promote_rule(::Type{GBMatrix}, ::Type{GBMatrixC}) = GBMatrix
41+
Base.promote_rule(::Type{GBMatrix}, ::Type{GBMatrixR}) = GBMatrix
42+
Base.promote_rule(::Type{GBMatrix}, ::Type{GBShallowMatrix}) = GBMatrix
43+
Base.promote_rule(::Type{GBMatrixC}, ::Type{GBMatrixC}) = GBMatrixC
44+
Base.promote_rule(::Type{GBMatrixR}, ::Type{GBMatrixR}) = GBMatrixR
45+
Base.promote_rule(::Type{GBMatrixC}, ::Type{GBMatrixR}) = GBMatrix
46+
Base.promote_rule(::Type{GBMatrixC}, ::Type{GBShallowMatrix}) = GBMatrixC
47+
Base.promote_rule(::Type{GBMatrixR}, ::Type{GBShallowMatrix}) = GBMatrixR
48+
49+
Base.promote_rule(::Type{G}, ::Type{<:AbstractGBVector}) where {G<:AbstractGBMatrix} = G
50+
Base.promote_rule(::Type{GBShallowMatrix}, ::Type{<:AbstractGBVector}) = GBMatrix
51+
52+
Base.promote_rule(::Type{<:AbstractGBVector}, ::Type{<:AbstractGBVector}) = GBVector
53+
Base.promote_rule(::Type{<:AbstractGBVector{T, F}}, ::Type{<:AbstractGBVector{T2, F2}}) where {T, F, T2, F2} =
54+
GBVector{promote_type(T, T2), promote_type(F, F2)}
55+
56+
function gbpromote_strip(A, B)
57+
if A isa Transpose{<:Any, <:AbstractGBVector} && B isa Transpose{<:Any, <:AbstractGBVector}
58+
return GBMatrix
59+
else
60+
return promote_type(strip_parameters(typeof(parent(A))), strip_parameters(typeof(parent(B))))
61+
end
62+
end
3263

3364
Base.IndexStyle(::AbstractGBArray) = IndexCartesian()
3465
Base.eltype(::AbstractGBArray{T, F}) where {T, F} = Union{T, F}

‎src/convert.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,7 @@ function Base.convert(::Type{M}, A::N; fill::F = getfill(A)) where {F, M<:Abstra
1616
newindices = _copytoraw.(indices)
1717
repack!()
1818
B = M(size(A, 1), size(A, 2); fill)
19-
unsafepack!(B, newindices..., newvalues, false; decrementindices = false)
19+
unsafepack!(B, newindices..., newvalues, false; decrementindices = false, order = storageorder(A))
2020
end
2121

2222
Base.convert(::Type{M}, A::M; fill = nothing) where {M<:AbstractGBArray} = A

‎src/gbmatrix.jl‎

Lines changed: 0 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -12,34 +12,3 @@ function Base.getindex(A::GBMatrixOrTranspose, v::AbstractVector)
1212
end
1313

1414
# Pack based constructors:
15-
function GBMatrix{T, F}(
16-
A::SparseVector;
17-
fill = defaultfill(F)
18-
) where {T, F}
19-
C = GBMatrix{T, F}(size(A, 1), 1; fill)
20-
return unsafepack!(C, _copytoraw(A)..., false)
21-
end
22-
GBMatrix{T}(
23-
A::SparseVector;
24-
fill::F = defaultfill(T)
25-
) where {T, F} = GBMatrix{T, F}(A; fill)
26-
GBMatrix(
27-
A::SparseVector{T};
28-
fill::F = defaultfill(T)
29-
) where {T, F} = GBMatrix{T, F}(A; fill)
30-
31-
function GBMatrix{T, F}(
32-
A::SparseMatrixCSC;
33-
fill = defaultfill(F)
34-
) where {T, F}
35-
C = GBMatrix{T, F}(size(A)...; fill)
36-
return unsafepack!(C, _copytoraw(A)..., false)
37-
end
38-
GBMatrix{T}(
39-
A::SparseMatrixCSC;
40-
fill::F = defaultfill(T)
41-
) where {T, F} = GBMatrix{T, F}(A; fill)
42-
GBMatrix(
43-
A::SparseMatrixCSC{T};
44-
fill::F = defaultfill(T)
45-
) where {T, F} = GBMatrix{T, F}(A; fill)

‎src/indexutils.jl‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -39,6 +39,9 @@ end
3939
# 1 <= length(szA | szB) <= 2
4040
# size checks should be done elsewhere.
4141
function _combinesizes(A, B)
42+
if A isa Transpose{<:Any, <:AbstractVector} && B isa Transpose{<:Any, <:AbstractVector}
43+
return size(A)
44+
end
4245
if (A isa AbstractVector && B isa AbstractMatrix) ||
4346
(B isa AbstractVector && A isa AbstractMatrix)
4447
return (size(A, 1), size(A, 2))

‎src/operations/broadcasts.jl‎

Lines changed: 9 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -48,14 +48,14 @@ function Base.similar(
4848
bc::Broadcast.Broadcasted{GBMatrixStyle},
4949
::Type{ElType}
5050
) where {ElType}
51-
return GBMatrix{ElType}(axes(bc))
51+
return GBMatrix{ElType}(axes(bc)) # Unfortunate default.
5252
end
5353

5454
function Base.similar(
5555
bc::Broadcast.Broadcasted{GBVectorStyle},
5656
::Type{ElType}
5757
) where {ElType}
58-
return GBVector{ElType}(axes(bc))
58+
return GBVector{ElType}(axes(bc)) # Okay default.
5959
end
6060

6161
#Find the modifying version of a function.
@@ -96,7 +96,13 @@ modifying(::typeof(emul)) = emul!
9696
add = defaultadd(f)
9797
return add(left, right, f)
9898
else
99-
return apply(f, left, right)
99+
leftscalar = !(left isa AbstractArray)
100+
rightscalar = !(right isa AbstractArray)
101+
if leftscalar || rightscalar
102+
return apply(f, left, right)
103+
end
104+
# TODO: WE NEED TO SUPPORT SARRAY ELTYPES HERE!!!
105+
return map(f, left, right)
100106
end
101107
end
102108
end

‎src/operations/concat.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -29,7 +29,7 @@ function Base.cat(tiles::VecOrMat{<:Union{AbstractGBMatrix{T, F}, AbstractGBVect
2929
end
3030
sz = (tiles isa AbstractArray && ncols == 1) ? (nrows,) : (nrows, ncols)
3131

32-
C = similar(tiles[1], t, sz; fill)
32+
C = similar(tiles[1], t, sz; fill) # TODO: FIXME, we want this to use promotion, but it's complicated.
3333
return cat!(C, tiles)
3434
end
3535

‎src/operations/ewise.jl‎

Lines changed: 31 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -77,9 +77,10 @@ function emul(
7777
accum = nothing,
7878
desc = nothing
7979
)
80-
t = inferbinarytype(parent(A), parent(B), op)
81-
82-
C = similar(A, t, _combinesizes(A, B); fill=_promotefill(parent(A), parent(B), op))
80+
T = inferbinarytype(parent(A), parent(B), op)
81+
fill=_promotefill(parent(A), parent(B), op)
82+
M = gbpromote_strip(A, B)
83+
C = M{T}(_combinesizes(A, B); fill)
8384
return emul!(C, A, B, op; mask, accum, desc)
8485
end
8586

@@ -163,8 +164,10 @@ function eadd(
163164
accum = nothing,
164165
desc = nothing
165166
)
166-
t = inferbinarytype(parent(A), parent(B), op)
167-
C = similar(A, t, _combinesizes(A, B); fill=_promotefill(parent(A), parent(B), op))
167+
T = inferbinarytype(parent(A), parent(B), op)
168+
fill=_promotefill(parent(A), parent(B), op)
169+
M = gbpromote_strip(A, B)
170+
C = M{T}(_combinesizes(A, B); fill)
168171
return eadd!(C, A, B, op; mask, accum, desc)
169172
end
170173

@@ -250,10 +253,18 @@ function eunion(
250253
desc = nothing
251254
) where {T, U}
252255
t = inferbinarytype(parent(A), parent(B), op)
253-
C = similar(A, t, _combinesizes(A, B); fill=_promotefill(parent(A), parent(B), op))
256+
fill=_promotefill(parent(A), parent(B), op)
257+
M = gbpromote_strip(A, B)
258+
C = M{t}(_combinesizes(A, B); fill)
254259
return eunion!(C, A, α, B, β, op; mask, accum, desc)
255260
end
256261

262+
eunion(
263+
A::GBArrayOrTranspose, α, B::GBArrayOrTranspose, β, op = +;
264+
mask = nothing, accum = nothing, desc = nothing
265+
) = eunion(A, convert(storedeltype(A), α), B, convert(storedeltype(B), β), op; mask, accum, desc)
266+
267+
257268
function Base.:+(A::GBArrayOrTranspose, B::GBArrayOrTranspose)
258269
eadd(A, B, +)
259270
end
@@ -299,4 +310,17 @@ eadd(A::VecMatOrTrans, B::GBArrayOrTranspose, op = +; kwargs...) =
299310
eadd(A::GBArrayOrTranspose, B::VecMatOrTrans, op = +; kwargs...) =
300311
@_densepack B eadd(A, B, op; kwargs...)
301312
eadd(A::VecMatOrTrans, B::VecMatOrTrans, op = +; kwargs...) =
302-
@_densepack A B eadd(A, B, op; kwargs...)
313+
@_densepack A B eadd(A, B, op; kwargs...)
314+
315+
eunion!(C::GBVecOrMat, A::VecMatOrTrans, α, B::GBArrayOrTranspose, β, op = +; kwargs...) =
316+
@_densepack A eunion!(C, A, α, B, β, op; kwargs...)
317+
eunion!(C::GBVecOrMat, A::GBArrayOrTranspose, α, B::VecMatOrTrans, β, op = +; kwargs...) =
318+
@_densepack B eunion!(C, A, α, B, β, op; kwargs...)
319+
eunion!(C::GBVecOrMat, A::VecMatOrTrans, α, B::VecMatOrTrans, β, op = +; kwargs...) =
320+
@_densepack A B eunion!(C, A, α, B, β, op; kwargs...)
321+
eunion(A::VecMatOrTrans, α, B::GBArrayOrTranspose, β, op = +; kwargs...) =
322+
@_densepack A eunion(A, α, B, β, op; kwargs...)
323+
eunion(A::GBArrayOrTranspose, α, B::VecMatOrTrans, β, op = +; kwargs...) =
324+
@_densepack B eunion(A, α, B, β, op; kwargs...)
325+
eunion(A::VecMatOrTrans, α, B::VecMatOrTrans, β, op = +; kwargs...) =
326+
@_densepack A B eunion(A, B, op; kwargs...)

‎src/operations/kronecker.jl‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -47,8 +47,10 @@ function LinearAlgebra.kron(
4747
accum = nothing,
4848
desc = nothing
4949
)
50-
t = inferbinarytype(parent(A), parent(B), op)
51-
C = similar(A, t, (size(A, 1) * size(B, 1), size(A, 2) * size(B, 2)); fill = _promotefill(parent(A), parent(B), op))
50+
T = inferbinarytype(parent(A), parent(B), op)
51+
fill=_promotefill(parent(A), parent(B), op)
52+
M = gbpromote_strip(A, B)
53+
C = M{T}(size(A, 1) * size(B, 1), size(A, 2) * size(B, 2); fill)
5254
kron!(C, A, B, op; mask, accum, desc)
5355
return C
5456
end

‎src/operations/mul.jl‎

Lines changed: 6 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -91,16 +91,17 @@ function Base.:*(
9191
accum = nothing,
9292
desc = nothing
9393
)
94-
t = inferbinarytype(parent(A), parent(B), op)
94+
T = inferbinarytype(parent(A), parent(B), op)
9595
fill = _promotefill(parent(A), parent(B), op)
9696
if A isa GBMatrixOrTranspose && B isa AbstractGBVector
97-
C = similar(A, t, size(A, 1); fill)
97+
C = similar(A, T, size(A, 1); fill)
9898
elseif A isa AbstractGBVector && B isa GBMatrixOrTranspose
99-
C = similar(A, t, size(B, 2); fill)
99+
C = similar(A, T, size(B, 2); fill)
100100
elseif A isa Transpose{<:Any, <:AbstractGBVector} && B isa AbstractGBVector
101-
C = similar(A, t, 1; fill)
101+
C = similar(A, T, 1; fill)
102102
else
103-
C = similar(A, t, (size(A, 1), size(B, 2)); fill)
103+
M = gbpromote_strip(A, B)
104+
C = M{T}((size(A, 1), size(B, 2)); fill)
104105
end
105106
mul!(C, A, B, op; mask, accum, desc)
106107
return C

‎src/operations/transpose.jl‎

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -66,6 +66,7 @@ end
6666

6767
#This is ok per the GraphBLAS Slack channel. Should change its effect on Complex input.
6868
LinearAlgebra.adjoint(A::GBVecOrMat) = transpose(A)
69+
LinearAlgebra.adjoint(A::GBVecOrMat{<:Complex}) = transpose(map(conj, A))
6970

7071
Base.unsafe_convert(::Type{Ptr{T}}, A::LinearAlgebra.AdjOrTrans{<:Any, <:AbstractGBArray}) where {T} =
7172
throw(ArgumentError("Cannot convert $(typeof(A)) directly to a pointer. Please use copy."))
@@ -96,4 +97,7 @@ Apply a mask to matrix `A`.
9697
"""
9798
function mask(A::GBArrayOrTranspose, mask; desc = nothing, replace_output = true)
9899
return mask!(similar(A), A, mask; desc)
99-
end
100+
end
101+
102+
# Gross, and should be better. But no real way around it currently:
103+

0 commit comments

Comments
 (0)