From cc6509482a12ab96b23d64a453740986706a6197 Mon Sep 17 00:00:00 2001 From: birkmichael <91692218+birkmichael@users.noreply.github.com> Date: Mon, 4 Oct 2021 18:33:09 +0300 Subject: [PATCH 01/21] calling similar on AbstractCuSparseArray with dims Calling similar() on a variable var of type : Date: Mon, 4 Oct 2021 22:18:14 +0300 Subject: [PATCH 02/21] robust handling of similar CSC/CSR Added more robust handling of applying similar on a CuSparseMatrixCSC/CSR, following the logic used in SparseArrays.jl --- lib/cusparse/array.jl | 63 ++++++++++++++++++++++++++++++++++++++----- 1 file changed, 57 insertions(+), 6 deletions(-) diff --git a/lib/cusparse/array.jl b/lib/cusparse/array.jl index 473cc7051b..9b928bb253 100644 --- a/lib/cusparse/array.jl +++ b/lib/cusparse/array.jl @@ -200,10 +200,59 @@ Base.similar(Mat::CuSparseMatrixCSC) = CuSparseMatrixCSC(copy(Mat.colPtr), copy( Base.similar(Mat::CuSparseMatrixCSR) = CuSparseMatrixCSR(copy(Mat.rowPtr), copy(Mat.colVal), similar(nonzeros(Mat)), Mat.dims) Base.similar(Mat::CuSparseMatrixBSR) = CuSparseMatrixBSR(copy(Mat.rowPtr), copy(Mat.colVal), similar(nonzeros(Mat)), Mat.blockDim, Mat.dir, nnz(Mat), Mat.dims) -Base.similar(a::AbstractCuSparseArray{Tv, <:Any, <:Any}, dims::Base.Dims{N}) where {Tv,N} = - CuArray{Tv,N}(undef, dims) -Base.similar(a::AbstractCuSparseArray{<:Any, <:Any, <:Any}, ::Type{T}, dims::Base.Dims{N}) where {T,N} = - CuArray{T,N}(undef, dims) +## similar for CSC,CSR with dims, eltype arguments, adapted from SparseArrays.jl. NOTE: calling similar() on a wrapped CuSparseMatrixBSR/COO will result in a CuArray. +# parent method for similar that preserves stored-entry structure (for when new and old dims match) +function _cusparsesimilar(S::CuSparseMatrixCSC, ::Type{TvNew}, ::Type{TiNew}) where {TvNew,TiNew} + newcolptr = copyto!(similar(getcolptr(S), TiNew), getcolptr(S)) + newrowval = copyto!(similar(rowvals(S), TiNew), rowvals(S)) + return CuSparseMatrixCSC(newcolptr, newrowval, similar(nonzeros(S), TvNew), size(S)) +end +function _cusparsesimilar(S::CuSparseMatrixCSR, ::Type{TvNew}, ::Type{TiNew}) where {TvNew,TiNew} + newrowptr = copyto!(similar(getrowptr(S), TiNew), getrowptr(S)) + newcolval = copyto!(similar(colvals(S), TiNew), colvals(S)) + return CuSparseMatrixCSR(newrowptr, newcolval, similar(nonzeros(S), TvNew), size(S)) +end +# parent methods for similar that preserves only storage space (for when new and old dims differ) +_cusparsesimilar(S::CuSparseMatrixCSC, ::Type{TvNew}, ::Type{TiNew}, dims::Dims{2}) where {TvNew,TiNew} = + CuSparseMatrixCSC(fill(one(TiNew), last(dims)+1), similar(rowvals(S), TiNew), similar(nonzeros(S), TvNew),dims) +_cusparsesimilar(S::CuSparseMatrixCSR, ::Type{TvNew}, ::Type{TiNew}, dims::Dims{2}) where {TvNew,TiNew} = + CuSparseMatrixCSR(fill(one(TiNew), last(dims)+1), similar(colvals(S), TiNew), similar(nonzeros(S), TvNew),dims) +# parent method for similar that allocates an empty sparse vector (when new dims are single) +_cusparsesimilar(S::CuSparseMatrixCSC, ::Type{TvNew}, ::Type{TiNew}, dims::Dims{1}) where {TvNew,TiNew} = + CuSparseVector(similar(rowvals(S), TiNew, 0), similar(nonzeros(S), TvNew, 0),dims) +_cusparsesimilar(S::CuSparseMatrixCSR, ::Type{TvNew}, ::Type{TiNew}, dims::Dims{1}) where {TvNew,TiNew} = + CuSparseVector(similar(colvals(S), TiNew, 0), similar(nonzeros(S), TvNew, 0),dims) + +""" +Utility union type of [`CuSparseMatrixCSC`](@ref), [`CuSparseMatrixCSR`](@ref), for which similar is implemented +""" +const CuSparseMatrixCSCR{Tv, Ti} = Union{ + CuSparseMatrixCSC{Tv, Ti}, + CuSparseMatrixCSR{Tv, Ti} +} +# The following methods hook into the AbstractArray similar hierarchy. The first method +# covers similar(A[, Tv]) calls, which preserve stored-entry structure, and the latter +# methods cover similar(A[, Tv], shape...) calls, which preserve storage space when the shape +# calls for a two-dimensional result. +similar(S::CuSparseMatrixCSCR{<:Any,Ti}, ::Type{TvNew}) where {Ti,TvNew} = _cusparsesimilar(S, TvNew, Ti) +similar(S::CuSparseMatrixCSCR{<:Any,Ti}, ::Type{TvNew}, dims::Union{Dims{1},Dims{2}}) where {Ti,TvNew} = + _cusparsesimilar(S, TvNew, Ti, dims) +# The following methods cover similar(A, Tv, Ti[, shape...]) calls, which specify the +# result's index type in addition to its entry type, and aren't covered by the hooks above. +# The calls without shape again preserve stored-entry structure, whereas those with shape +# preserve storage space when the shape calls for a two-dimensional result. +Base.similar(S::CuSparseMatrixCSCR, ::Type{TvNew}, ::Type{TiNew}) where{TvNew,TiNew} = + _cusparsesimilar(S, TvNew, TiNew) +Base.similar(S::CuSparseMatrixCSCR, ::Type{TvNew}, ::Type{TiNew}, dims::Union{Dims{1},Dims{2}}) where {TvNew,TiNew} = + _cusparsesimilar(S, TvNew, TiNew, dims) +Base.similar(S::CuSparseMatrixCSCR, ::Type{TvNew}, ::Type{TiNew}, m::Integer) where {TvNew,TiNew} = + _cusparsesimilar(S, TvNew, TiNew, (m,)) +Base.similar(S::CuSparseMatrixCSCR, ::Type{TvNew}, ::Type{TiNew}, m::Integer, n::Integer) where {TvNew,TiNew} = + _cusparsesimilar(S, TvNew, TiNew, (m, n)) + +# Handles cases where CSC/CSR are reshaped to more than two dimensions, and all cases above for COO/BSR. +Base.similar(a::AbstractCuSparseArray, ::Type{T}, dims::Base.Dims{N}) where {T,N} = +CuArray{T,N}(undef, dims) ## array interface @@ -242,9 +291,11 @@ Base.eltype(g::CuSparseMatrix{T}) where T = T SparseArrays.nnz(g::AbstractCuSparseArray) = g.nnz SparseArrays.nonzeros(g::AbstractCuSparseArray) = g.nzVal - SparseArrays.nonzeroinds(g::AbstractCuSparseVector) = g.iPtr - + +SparseArrays.getcolptr(g::CuSparseMatrixCSC) = g.colPtr +SparseArrays.getrowptr(g::CuSparseMatrixCSR) = g.rowPtr +SparseArrays.colvals(g::CuSparseMatrixCSR) = g.colVal SparseArrays.rowvals(g::CuSparseMatrixCSC) = g.rowVal LinearAlgebra.issymmetric(M::Union{CuSparseMatrixCSC,CuSparseMatrixCSR}) = false From 53724cea0c2f018ef2eb1568b3a0de991dae0372 Mon Sep 17 00:00:00 2001 From: birkmichael <91692218+birkmichael@users.noreply.github.com> Date: Mon, 4 Oct 2021 22:43:58 +0300 Subject: [PATCH 03/21] fixed getrowptr,colvals two methods CSR were originally added as extensions to SparseArrays.jl, but SparseArrays.jl does not support CSR. --- lib/cusparse/array.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/lib/cusparse/array.jl b/lib/cusparse/array.jl index 9b928bb253..3ae148ab85 100644 --- a/lib/cusparse/array.jl +++ b/lib/cusparse/array.jl @@ -294,9 +294,10 @@ SparseArrays.nonzeros(g::AbstractCuSparseArray) = g.nzVal SparseArrays.nonzeroinds(g::AbstractCuSparseVector) = g.iPtr SparseArrays.getcolptr(g::CuSparseMatrixCSC) = g.colPtr -SparseArrays.getrowptr(g::CuSparseMatrixCSR) = g.rowPtr -SparseArrays.colvals(g::CuSparseMatrixCSR) = g.colVal SparseArrays.rowvals(g::CuSparseMatrixCSC) = g.rowVal +getrowptr(g::CuSparseMatrixCSR) = g.rowPtr +colvals(g::CuSparseMatrixCSR) = g.colVal + LinearAlgebra.issymmetric(M::Union{CuSparseMatrixCSC,CuSparseMatrixCSR}) = false LinearAlgebra.ishermitian(M::Union{CuSparseMatrixCSC,CuSparseMatrixCSR}) = false From 6a8b64bd682ab2b56f045bd05df598fce0fe1e02 Mon Sep 17 00:00:00 2001 From: birkmichael <91692218+birkmichael@users.noreply.github.com> Date: Mon, 4 Oct 2021 22:59:24 +0300 Subject: [PATCH 04/21] forgot Base. prefix to several instances of similar --- lib/cusparse/array.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/cusparse/array.jl b/lib/cusparse/array.jl index 3ae148ab85..d14c89fba2 100644 --- a/lib/cusparse/array.jl +++ b/lib/cusparse/array.jl @@ -234,8 +234,8 @@ const CuSparseMatrixCSCR{Tv, Ti} = Union{ # covers similar(A[, Tv]) calls, which preserve stored-entry structure, and the latter # methods cover similar(A[, Tv], shape...) calls, which preserve storage space when the shape # calls for a two-dimensional result. -similar(S::CuSparseMatrixCSCR{<:Any,Ti}, ::Type{TvNew}) where {Ti,TvNew} = _cusparsesimilar(S, TvNew, Ti) -similar(S::CuSparseMatrixCSCR{<:Any,Ti}, ::Type{TvNew}, dims::Union{Dims{1},Dims{2}}) where {Ti,TvNew} = +Base.similar(S::CuSparseMatrixCSCR{<:Any,Ti}, ::Type{TvNew}) where {Ti,TvNew} = _cusparsesimilar(S, TvNew, Ti) +Base.similar(S::CuSparseMatrixCSCR{<:Any,Ti}, ::Type{TvNew}, dims::Union{Dims{1},Dims{2}}) where {Ti,TvNew} = _cusparsesimilar(S, TvNew, Ti, dims) # The following methods cover similar(A, Tv, Ti[, shape...]) calls, which specify the # result's index type in addition to its entry type, and aren't covered by the hooks above. From 2620166fb5db97739775aa9ab15b744bb60fed46 Mon Sep 17 00:00:00 2001 From: birkmichael <91692218+birkmichael@users.noreply.github.com> Date: Mon, 4 Oct 2021 23:24:47 +0300 Subject: [PATCH 05/21] used fill instead of CUDA.fill in matrix initialization --- lib/cusparse/array.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/cusparse/array.jl b/lib/cusparse/array.jl index d14c89fba2..b932f9afa1 100644 --- a/lib/cusparse/array.jl +++ b/lib/cusparse/array.jl @@ -214,9 +214,9 @@ function _cusparsesimilar(S::CuSparseMatrixCSR, ::Type{TvNew}, ::Type{TiNew}) wh end # parent methods for similar that preserves only storage space (for when new and old dims differ) _cusparsesimilar(S::CuSparseMatrixCSC, ::Type{TvNew}, ::Type{TiNew}, dims::Dims{2}) where {TvNew,TiNew} = - CuSparseMatrixCSC(fill(one(TiNew), last(dims)+1), similar(rowvals(S), TiNew), similar(nonzeros(S), TvNew),dims) + CuSparseMatrixCSC(CUDA.fill(one(TiNew), last(dims)+1), similar(rowvals(S), TiNew), similar(nonzeros(S), TvNew),dims) _cusparsesimilar(S::CuSparseMatrixCSR, ::Type{TvNew}, ::Type{TiNew}, dims::Dims{2}) where {TvNew,TiNew} = - CuSparseMatrixCSR(fill(one(TiNew), last(dims)+1), similar(colvals(S), TiNew), similar(nonzeros(S), TvNew),dims) + CuSparseMatrixCSR(CUDA.fill(one(TiNew), last(dims)+1), similar(colvals(S), TiNew), similar(nonzeros(S), TvNew),dims) # parent method for similar that allocates an empty sparse vector (when new dims are single) _cusparsesimilar(S::CuSparseMatrixCSC, ::Type{TvNew}, ::Type{TiNew}, dims::Dims{1}) where {TvNew,TiNew} = CuSparseVector(similar(rowvals(S), TiNew, 0), similar(nonzeros(S), TvNew, 0),dims) From 05c85ef32061461b4767c62d01035c3e13a7d0b5 Mon Sep 17 00:00:00 2001 From: birkmichael <91692218+birkmichael@users.noreply.github.com> Date: Tue, 5 Oct 2021 21:54:06 +0300 Subject: [PATCH 06/21] tests for similar, unfinished --- test/cusparse/array.jl | 208 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 208 insertions(+) create mode 100644 test/cusparse/array.jl diff --git a/test/cusparse/array.jl b/test/cusparse/array.jl new file mode 100644 index 0000000000..6ca8658377 --- /dev/null +++ b/test/cusparse/array.jl @@ -0,0 +1,208 @@ +using CUDA.CUSPARSE + +@testset "similar" begin + @testset "$f(A)±$h(B) $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64], + f in (identity, transpose), #adjoint), + (h,dimPtr,dimVal) in ((CuSparseMatrixCSR,:rowPtr,:colVal), (CuSparseMatrixCSC,:colPtr,:rowVal)) + n = 10 + A = sprand(elty, n, n, rand()) + dA = f(A) + + C_simple = similar(dA) + C_dims = similar(dA,(n,n+1)) + C_eltype = similar(dA,elty) + C_full = similar(dA,elty,(n,n+1)) + + @test typeof(C_simple) ≈ typeof(dA) + @test getproperty(C_simple,dimPtr) ≈ getproperty(dA,dimPtr) #similar without arguments perserves structure + @test getproperty(C_simple,dimVal) ≈ getproperty(dA,dimVal) + @test getproperty(C_simple,dims) ≈ getproperty(dA,dims) + + @test typeof(C_simple) ≈ typeof(dA) + @test getproperty(C_simple,dimPtr) ≈ getproperty(dA,dimPtr) #similar without arguments perserves structure + @test getproperty(C_simple,dimVal) ≈ getproperty(dA,dimVal) + @test getproperty(C_dims,dims) ≈ (n,n+1) + for propertyname in propertynames(dA) + @test length(getproperty(C_simple,z)) ≈ length(getproperty(dA,z)) + end + @test C ≈ collect(dC) + + C = f(A) - h(B) + dC = f(dA) - h(dB) + @test C ≈ collect(dC) + end + + @testset "A±$f(B) $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64], + f in (CuSparseMatrixCSR, CuSparseMatrixCSC, x->CuSparseMatrixBSR(x,1)) + n = 10 + A = sprand(elty, n, n, rand()) + B = sprand(elty, n, n, rand()) + + dA = CuSparseMatrixCSR(A) + dB = CuSparseMatrixCSR(B) + + C = A + B + dC = dA + f(dB) + @test C ≈ collect(dC) + + C = B + A + dC = f(dB) + dA + @test C ≈ collect(dC) + + C = A - B + dC = dA - f(dB) + @test C ≈ collect(dC) + + C = B - A + dC = f(dB) - dA + @test C ≈ collect(dC) + end + + @testset "dense(A)$(op)sparse(B) $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64], + op in [+, -] + n = 10 + A = rand(elty, n, n) + B = sprand(elty, n, n, rand()) + + dA = CuArray(A) + dB = CuSparseMatrixCSR(B) + + C = op(A, B) + dC = op(dA, dB) + @test C ≈ collect(dC) + @test dC isa CuMatrix{elty} + + C = op(B, A) + dC = op(dB, dA) + @test C ≈ collect(dC) + @test dC isa CuMatrix{elty} + end + + @testset "$f(A)*b $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64], + f in (identity, transpose, adjoint) + n = 10 + alpha = rand() + beta = rand() + A = sprand(elty, n, n, rand()) + b = rand(elty, n) + c = rand(elty, n) + alpha = beta = 1.0 + c = zeros(elty, n) + + dA = CuSparseMatrixCSR(A) + db = CuArray(b) + dc = CuArray(c) + + # test with empty inputs + @test Array(dA * CUDA.zeros(elty, n, 0)) == zeros(elty, n, 0) + + mul!(c, f(A), b, alpha, beta) + mul!(dc, f(dA), db, alpha, beta) + @test c ≈ collect(dc) + + A = A + transpose(A) + dA = CuSparseMatrixCSR(A) + + @assert issymmetric(A) + mul!(c, f(Symmetric(A)), b, alpha, beta) + mul!(dc, f(Symmetric(dA)), db, alpha, beta) + @test c ≈ collect(dc) + end + + @testset "$f(A)*$h(B) $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64], + f in (identity, transpose, adjoint), + h in (identity, transpose, adjoint) + if CUSPARSE.version() <= v"10.3.1" && + (h ∈ (transpose, adjoint) && f != identity) || + (h == adjoint && elty <: Complex) + # These combinations are not implemented in CUDA10 + continue + end + + n = 10 + alpha = rand() + beta = rand() + A = sprand(elty, n, n, rand()) + B = rand(elty, n, n) + C = rand(elty, n, n) + + dA = CuSparseMatrixCSR(A) + dB = CuArray(B) + dC = CuArray(C) + + mul!(C, f(A), h(B), alpha, beta) + mul!(dC, f(dA), h(dB), alpha, beta) + @test C ≈ collect(dC) + + A = A + transpose(A) + dA = CuSparseMatrixCSR(A) + + @assert issymmetric(A) + mul!(C, f(Symmetric(A)), h(B), alpha, beta) + mul!(dC, f(Symmetric(dA)), h(dB), alpha, beta) + @test C ≈ collect(dC) + end + + @testset "dense(A)*sparse(B) $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64] + n = 10 + A = rand(elty, n, n) + B = sprand(elty, n, n, rand()) + + dA = CuArray(A) + dB = CUSPARSE.CuSparseMatrixCSR(B) + + C = A * B + dC = dA * dB + @test C ≈ collect(dC) + @test dC isa CuMatrix{elty} + + C = B * A + dC = dB * dA + @test C ≈ collect(dC) + @test dC isa CuMatrix{elty} + + C = B * B + dC = dB * dB + @test C ≈ collect(dC) + @test dC isa CuMatrix{elty} + end + + @testset "issue #1095 ($elty)" for elty in [Float32, Float64, ComplexF32, ComplexF64] + # Test non-square matrices + n, m, p = 10, 20, 4 + alpha = rand() + beta = rand() + A = sprand(elty, n, m, rand()) + B = rand(elty, m, p) + C = rand(elty, n, p) + + dA = CuSparseMatrixCSR(A) + dB = CuArray(B) + dC = CuArray(C) + + mul!(C, A, B, alpha, beta) + mul!(dC, dA, dB, alpha, beta) + @test C ≈ collect(dC) + + mul!(B, transpose(A), C, alpha, beta) + mul!(dB, transpose(dA), dC, alpha, beta) + @test B ≈ collect(dB) + + mul!(B, adjoint(A), C, alpha, beta) + mul!(dB, adjoint(dA), dC, alpha, beta) + @test B ≈ collect(dB) + end + + @testset "CuSparseMatrixCSR($f) $elty" for f in [transpose, adjoint], elty in [Float32, ComplexF32] + S = f(sprand(elty, 10, 10, 0.1)) + @test SparseMatrixCSC(CuSparseMatrixCSR(S)) ≈ S + + S = sprand(elty, 10, 10, 0.1) + T = f(CuSparseMatrixCSR(S)) + @test SparseMatrixCSC(CuSparseMatrixCSC(T)) ≈ f(S) + + S = sprand(elty, 10, 10, 0.1) + T = f(CuSparseMatrixCSC(S)) + @test SparseMatrixCSC(CuSparseMatrixCSR(T)) ≈ f(S) + end +end From 309f02e05916c62be407083dbd735ab8b6420e90 Mon Sep 17 00:00:00 2001 From: birkmichael <91692218+birkmichael@users.noreply.github.com> Date: Wed, 6 Oct 2021 15:27:20 +0300 Subject: [PATCH 07/21] Created test set for 'similar' --- test/cusparse/array.jl | 289 +++++++++++++++-------------------------- 1 file changed, 103 insertions(+), 186 deletions(-) diff --git a/test/cusparse/array.jl b/test/cusparse/array.jl index 6ca8658377..178df4064b 100644 --- a/test/cusparse/array.jl +++ b/test/cusparse/array.jl @@ -1,208 +1,125 @@ using CUDA.CUSPARSE - +using SparseArrays +using Test +typeSet = [Float32, Float64, ComplexF32, ComplexF64] +n = 10 @testset "similar" begin - @testset "$f(A)±$h(B) $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64], - f in (identity, transpose), #adjoint), + @testset "similar(A::$h{$elty},$newElty)" for elty in typeSet, + newElty in typeSet, (h,dimPtr,dimVal) in ((CuSparseMatrixCSR,:rowPtr,:colVal), (CuSparseMatrixCSC,:colPtr,:rowVal)) - n = 10 + A = sprand(elty, n, n, rand()) - dA = f(A) + dA = h(A) C_simple = similar(dA) C_dims = similar(dA,(n,n+1)) - C_eltype = similar(dA,elty) - C_full = similar(dA,elty,(n,n+1)) - - @test typeof(C_simple) ≈ typeof(dA) - @test getproperty(C_simple,dimPtr) ≈ getproperty(dA,dimPtr) #similar without arguments perserves structure - @test getproperty(C_simple,dimVal) ≈ getproperty(dA,dimVal) - @test getproperty(C_simple,dims) ≈ getproperty(dA,dims) - - @test typeof(C_simple) ≈ typeof(dA) - @test getproperty(C_simple,dimPtr) ≈ getproperty(dA,dimPtr) #similar without arguments perserves structure - @test getproperty(C_simple,dimVal) ≈ getproperty(dA,dimVal) - @test getproperty(C_dims,dims) ≈ (n,n+1) - for propertyname in propertynames(dA) - @test length(getproperty(C_simple,z)) ≈ length(getproperty(dA,z)) - end - @test C ≈ collect(dC) - - C = f(A) - h(B) - dC = f(dA) - h(dB) - @test C ≈ collect(dC) - end - - @testset "A±$f(B) $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64], - f in (CuSparseMatrixCSR, CuSparseMatrixCSC, x->CuSparseMatrixBSR(x,1)) - n = 10 - A = sprand(elty, n, n, rand()) - B = sprand(elty, n, n, rand()) - - dA = CuSparseMatrixCSR(A) - dB = CuSparseMatrixCSR(B) - - C = A + B - dC = dA + f(dB) - @test C ≈ collect(dC) + C_eltype = similar(dA,newElty) + C_full = similar(dA,newElty,(n,n+1)) + @test typeof(C_simple) ≈ typeof(dA) + @test typeof(C_dims) ≈ typeof(dA) + @test typeof(C_eltype) ≈ typeof(dA) + @test typeof(C_full) ≈ typeof(dA) + + properties = Set([propertynames(dA)]); + + conserved_simple = Set([dimPtr,dimVal,dims,nnz]) + structure_conserved_simple = setdiff(properties,conserved_simple); + for propertyname in conserved_simple + @test getproperty(C_simple,propertyname) ≈ getproperty(dA,propertyname) + end + for propertyname in structure_conserved_simple + @test length(getproperty(C_simple,propertyname)) ≈ length(getproperty(dA,propertyname)) + @test eltype(getproperty(C_simple,propertyname)) ≈ eltype(getproperty(dA,propertyname)) + end - C = B + A - dC = f(dB) + dA - @test C ≈ collect(dC) + conserved_dims = Set([nnz]) + structure_conserved_dims = setdiff(properties,union(conserved_dims,dimVal,:dims)) + @test getproperty(C_dims,:dims) ≈ (n,n+1) + for propertyname in conserved_dims + @test getproperty(C_dims,propertyname) ≈ getproperty(dA,propertyname) + end + for propertyname in structure_conserved_dims + @test length(getproperty(C_dims,propertyname)) ≈ length(getproperty(dA,propertyname)) + @test eltype(getproperty(C_dims,propertyname)) ≈ eltype(getproperty(dA,propertyname)) + end - C = A - B - dC = dA - f(dB) - @test C ≈ collect(dC) + conserved_eltype = Set([nnz,dims,dimPtr,dimVal]) + structure_conserved_eltype = setdiff(properties,union(conserved_eltype,:nzVal)) + @test eltype(getproperty(C_eltype,:nzVal)) ≈ newElty + @test length(getproperty(C_eltype,:nzVal)) ≈ length(getproperty(dA,:nzVal)) + for propertyname in conserved_eltype + @test getproperty(C_dims,propertyname) ≈ getproperty(dA,propertyname) + end + for propertyname in structure_conserved_eltype + @test length(getproperty(C_dims,propertyname)) ≈ length(getproperty(dA,propertyname)) + @test eltype(getproperty(C_dims,propertyname)) ≈ eltype(getproperty(dA,propertyname)) + end - C = B - A - dC = f(dB) - dA - @test C ≈ collect(dC) + @test eltype(getproperty(C_full,:nzVal)) ≈ newElty + @test length(getproperty(C_full,:nzVal)) ≈ length(getproperty(dA,:nzVal)) + @test length(getproperty(C_full,:dimPtr)) ≈ length(getproperty(dA,:dimPtr)) + @test getproperty(C_dims,:nnz) ≈ getproperty(dA,:nnz) + @test getproperty(C_full,:dims) ≈ (n,n+1) end - - @testset "dense(A)$(op)sparse(B) $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64], - op in [+, -] - n = 10 - A = rand(elty, n, n) - B = sprand(elty, n, n, rand()) - - dA = CuArray(A) - dB = CuSparseMatrixCSR(B) - - C = op(A, B) - dC = op(dA, dB) - @test C ≈ collect(dC) - @test dC isa CuMatrix{elty} - C = op(B, A) - dC = op(dB, dA) - @test C ≈ collect(dC) - @test dC isa CuMatrix{elty} - end + @testset "similar(A::$f($h{$elty}),$newElty)" for elty in typeSet, + newElty in typeSet, + f in (transpose, x->reshape(x,n,n)), + (h,dimPtr,dimVal) in ((CuSparseMatrixCSR,:rowPtr,:colVal), (CuSparseMatrixCSC,:colPtr,:rowVal)) - @testset "$f(A)*b $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64], - f in (identity, transpose, adjoint) - n = 10 - alpha = rand() - beta = rand() A = sprand(elty, n, n, rand()) - b = rand(elty, n) - c = rand(elty, n) - alpha = beta = 1.0 - c = zeros(elty, n) - - dA = CuSparseMatrixCSR(A) - db = CuArray(b) - dc = CuArray(c) - - # test with empty inputs - @test Array(dA * CUDA.zeros(elty, n, 0)) == zeros(elty, n, 0) - - mul!(c, f(A), b, alpha, beta) - mul!(dc, f(dA), db, alpha, beta) - @test c ≈ collect(dc) - - A = A + transpose(A) - dA = CuSparseMatrixCSR(A) - - @assert issymmetric(A) - mul!(c, f(Symmetric(A)), b, alpha, beta) - mul!(dc, f(Symmetric(dA)), db, alpha, beta) - @test c ≈ collect(dc) - end - - @testset "$f(A)*$h(B) $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64], - f in (identity, transpose, adjoint), - h in (identity, transpose, adjoint) - if CUSPARSE.version() <= v"10.3.1" && - (h ∈ (transpose, adjoint) && f != identity) || - (h == adjoint && elty <: Complex) - # These combinations are not implemented in CUDA10 - continue + dA = f(h(A)) + + C_simple = similar(dA) + C_dims = similar(dA,(n,n+1)) + C_eltype = similar(dA,newElty) + C_full = similar(dA,newElty,(n,n+1)) + @test typeof(C_simple) ≈ typeof(dA) + @test typeof(C_dims) ≈ typeof(dA) + @test typeof(C_eltype) ≈ typeof(dA) + @test typeof(C_full) ≈ typeof(dA) + + properties = Set([propertynames(dA)]); + + conserved_simple = Set([nnz]) + structure_conserved_simple = setdiff(properties,conserved_simple); + for propertyname in conserved_simple + @test getproperty(C_simple,propertyname) ≈ getproperty(dA,propertyname) + end + for propertyname in structure_conserved_simple + @test length(getproperty(C_simple,propertyname)) ≈ length(getproperty(dA,propertyname)) + @test eltype(getproperty(C_simple,propertyname)) ≈ eltype(getproperty(dA,propertyname)) end - n = 10 - alpha = rand() - beta = rand() - A = sprand(elty, n, n, rand()) - B = rand(elty, n, n) - C = rand(elty, n, n) - - dA = CuSparseMatrixCSR(A) - dB = CuArray(B) - dC = CuArray(C) - - mul!(C, f(A), h(B), alpha, beta) - mul!(dC, f(dA), h(dB), alpha, beta) - @test C ≈ collect(dC) + conserved_dims = Set([nnz]) + structure_conserved_dims = setdiff(properties,union(conserved_dims,dimVal,:dims)) + @test getproperty(C_dims,:dims) ≈ (n,n+1) + for propertyname in conserved_dims + @test getproperty(C_dims,propertyname) ≈ getproperty(dA,propertyname) + end + for propertyname in structure_conserved_dims + @test length(getproperty(C_dims,propertyname)) ≈ length(getproperty(dA,propertyname)) + @test eltype(getproperty(C_dims,propertyname)) ≈ eltype(getproperty(dA,propertyname)) + end - A = A + transpose(A) - dA = CuSparseMatrixCSR(A) + conserved_eltype = Set([nnz,dims,dimPtr,dimVal]) + structure_conserved_eltype = setdiff(properties,union(conserved_eltype,:nzVal)) + @test eltype(getproperty(C_eltype,:nzVal)) ≈ newElty + @test length(getproperty(C_eltype,:nzVal)) ≈ length(getproperty(dA,:nzVal)) + for propertyname in conserved_eltype + @test getproperty(C_dims,propertyname) ≈ getproperty(dA,propertyname) + end + for propertyname in structure_conserved_eltype + @test length(getproperty(C_dims,propertyname)) ≈ length(getproperty(dA,propertyname)) + @test eltype(getproperty(C_dims,propertyname)) ≈ eltype(getproperty(dA,propertyname)) + end - @assert issymmetric(A) - mul!(C, f(Symmetric(A)), h(B), alpha, beta) - mul!(dC, f(Symmetric(dA)), h(dB), alpha, beta) - @test C ≈ collect(dC) + @test eltype(getproperty(C_full,:nzVal)) ≈ newElty + @test length(getproperty(C_full,:nzVal)) ≈ length(getproperty(dA,:nzVal)) + @test length(getproperty(C_full,:dimPtr)) ≈ length(getproperty(dA,:dimPtr)) + @test getproperty(C_dims,:nnz) ≈ getproperty(dA,:nnz) + @test getproperty(C_full,:dims) ≈ (n,n+1) end - @testset "dense(A)*sparse(B) $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64] - n = 10 - A = rand(elty, n, n) - B = sprand(elty, n, n, rand()) - - dA = CuArray(A) - dB = CUSPARSE.CuSparseMatrixCSR(B) - C = A * B - dC = dA * dB - @test C ≈ collect(dC) - @test dC isa CuMatrix{elty} - - C = B * A - dC = dB * dA - @test C ≈ collect(dC) - @test dC isa CuMatrix{elty} - - C = B * B - dC = dB * dB - @test C ≈ collect(dC) - @test dC isa CuMatrix{elty} - end - - @testset "issue #1095 ($elty)" for elty in [Float32, Float64, ComplexF32, ComplexF64] - # Test non-square matrices - n, m, p = 10, 20, 4 - alpha = rand() - beta = rand() - A = sprand(elty, n, m, rand()) - B = rand(elty, m, p) - C = rand(elty, n, p) - - dA = CuSparseMatrixCSR(A) - dB = CuArray(B) - dC = CuArray(C) - - mul!(C, A, B, alpha, beta) - mul!(dC, dA, dB, alpha, beta) - @test C ≈ collect(dC) - - mul!(B, transpose(A), C, alpha, beta) - mul!(dB, transpose(dA), dC, alpha, beta) - @test B ≈ collect(dB) - - mul!(B, adjoint(A), C, alpha, beta) - mul!(dB, adjoint(dA), dC, alpha, beta) - @test B ≈ collect(dB) - end - - @testset "CuSparseMatrixCSR($f) $elty" for f in [transpose, adjoint], elty in [Float32, ComplexF32] - S = f(sprand(elty, 10, 10, 0.1)) - @test SparseMatrixCSC(CuSparseMatrixCSR(S)) ≈ S - - S = sprand(elty, 10, 10, 0.1) - T = f(CuSparseMatrixCSR(S)) - @test SparseMatrixCSC(CuSparseMatrixCSC(T)) ≈ f(S) - - S = sprand(elty, 10, 10, 0.1) - T = f(CuSparseMatrixCSC(S)) - @test SparseMatrixCSC(CuSparseMatrixCSR(T)) ≈ f(S) - end end From 6ec07de477237e6b3fa5a46611f152fda3ecf805 Mon Sep 17 00:00:00 2001 From: birkmichael <91692218+birkmichael@users.noreply.github.com> Date: Wed, 6 Oct 2021 15:29:09 +0300 Subject: [PATCH 08/21] imported some SparseArrays methods --- lib/cusparse/array.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/cusparse/array.jl b/lib/cusparse/array.jl index b932f9afa1..9739ab9109 100644 --- a/lib/cusparse/array.jl +++ b/lib/cusparse/array.jl @@ -6,7 +6,7 @@ export CuSparseMatrixCSC, CuSparseMatrixCSR, CuSparseMatrixBSR, CuSparseMatrixCO CuSparseVector using LinearAlgebra: BlasFloat -using SparseArrays: nonzeroinds, dimlub +using SparseArrays: nonzeroinds, dimlub, getcolptr, rowvals abstract type AbstractCuSparseArray{Tv, Ti, N} <: AbstractSparseArray{Tv, Ti, N} end const AbstractCuSparseVector{Tv, Ti} = AbstractCuSparseArray{Tv, Ti, 1} From fb4b2535189bde26ec9048942e6f9066f4bc7898 Mon Sep 17 00:00:00 2001 From: birkmichael <91692218+birkmichael@users.noreply.github.com> Date: Wed, 6 Oct 2021 18:26:19 +0300 Subject: [PATCH 09/21] Fixed CSR rowPtr setting in 'similar' --- lib/cusparse/array.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/cusparse/array.jl b/lib/cusparse/array.jl index 9739ab9109..711a178fad 100644 --- a/lib/cusparse/array.jl +++ b/lib/cusparse/array.jl @@ -216,7 +216,7 @@ end _cusparsesimilar(S::CuSparseMatrixCSC, ::Type{TvNew}, ::Type{TiNew}, dims::Dims{2}) where {TvNew,TiNew} = CuSparseMatrixCSC(CUDA.fill(one(TiNew), last(dims)+1), similar(rowvals(S), TiNew), similar(nonzeros(S), TvNew),dims) _cusparsesimilar(S::CuSparseMatrixCSR, ::Type{TvNew}, ::Type{TiNew}, dims::Dims{2}) where {TvNew,TiNew} = - CuSparseMatrixCSR(CUDA.fill(one(TiNew), last(dims)+1), similar(colvals(S), TiNew), similar(nonzeros(S), TvNew),dims) + CuSparseMatrixCSR(CUDA.fill(one(TiNew), first(dims)+1), similar(colvals(S), TiNew), similar(nonzeros(S), TvNew),dims) # parent method for similar that allocates an empty sparse vector (when new dims are single) _cusparsesimilar(S::CuSparseMatrixCSC, ::Type{TvNew}, ::Type{TiNew}, dims::Dims{1}) where {TvNew,TiNew} = CuSparseVector(similar(rowvals(S), TiNew, 0), similar(nonzeros(S), TvNew, 0),dims) From df04acb045dddbb76335754625d336453954f6c6 Mon Sep 17 00:00:00 2001 From: birkmichael <91692218+birkmichael@users.noreply.github.com> Date: Wed, 6 Oct 2021 19:17:16 +0300 Subject: [PATCH 10/21] tests for 'similar' for CSC,CSR unclear what is supported in CUDA.jl for BSR,COO, so they remain untested. --- test/cusparse/array.jl | 154 +++++++++++++++++++++++------------------ 1 file changed, 86 insertions(+), 68 deletions(-) diff --git a/test/cusparse/array.jl b/test/cusparse/array.jl index 178df4064b..30b23a84a3 100644 --- a/test/cusparse/array.jl +++ b/test/cusparse/array.jl @@ -4,7 +4,7 @@ using Test typeSet = [Float32, Float64, ComplexF32, ComplexF64] n = 10 @testset "similar" begin - @testset "similar(A::$h{$elty},$newElty)" for elty in typeSet, + @testset "similar(A::$h{$elty},Tv::$newElty)" for elty in typeSet, newElty in typeSet, (h,dimPtr,dimVal) in ((CuSparseMatrixCSR,:rowPtr,:colVal), (CuSparseMatrixCSC,:colPtr,:rowVal)) @@ -15,51 +15,62 @@ n = 10 C_dims = similar(dA,(n,n+1)) C_eltype = similar(dA,newElty) C_full = similar(dA,newElty,(n,n+1)) - @test typeof(C_simple) ≈ typeof(dA) - @test typeof(C_dims) ≈ typeof(dA) - @test typeof(C_eltype) ≈ typeof(dA) - @test typeof(C_full) ≈ typeof(dA) + @test typeof(C_simple) == typeof(dA) + @test typeof(C_dims) == typeof(dA) + @test (typeof(C_eltype) == typeof(dA) && elty==newElty) || ((typeof(C_eltype) != typeof(dA) && elty!=newElty)) + @test (typeof(C_full) == typeof(dA) && elty==newElty) || ((typeof(C_full) != typeof(dA) && elty!=newElty)) - properties = Set([propertynames(dA)]); + properties = Set(propertynames(dA)); - conserved_simple = Set([dimPtr,dimVal,dims,nnz]) + conserved_simple = Set([dimPtr,dimVal,:dims,:nnz]) structure_conserved_simple = setdiff(properties,conserved_simple); - for propertyname in conserved_simple - @test getproperty(C_simple,propertyname) ≈ getproperty(dA,propertyname) + @testset "similar(A::$h{$elty},Tv::$newElty) $propertyname simple conserved" for propertyname in conserved_simple + @test getproperty(C_simple,propertyname) == getproperty(dA,propertyname) end - for propertyname in structure_conserved_simple - @test length(getproperty(C_simple,propertyname)) ≈ length(getproperty(dA,propertyname)) - @test eltype(getproperty(C_simple,propertyname)) ≈ eltype(getproperty(dA,propertyname)) + + @testset "similar(A::$h{$elty},Tv::$newElty) $propertyname simple structure conserved" for propertyname in structure_conserved_simple + @test length(getproperty(C_simple,propertyname)) == length(getproperty(dA,propertyname)) + @test eltype(getproperty(C_simple,propertyname)) == eltype(getproperty(dA,propertyname)) end - conserved_dims = Set([nnz]) - structure_conserved_dims = setdiff(properties,union(conserved_dims,dimVal,:dims)) - @test getproperty(C_dims,:dims) ≈ (n,n+1) - for propertyname in conserved_dims - @test getproperty(C_dims,propertyname) ≈ getproperty(dA,propertyname) + conserved_dims = Set([:nnz]) + if h==CuSparseMatrixCSR # Making the array one column longer increases colPtr length but not rowPtr length + structure_conserved_dims = setdiff(properties,union(conserved_dims,Set([dimVal,:dims]))) + else #CSC + structure_conserved_dims = setdiff(properties,union(conserved_dims,Set([dimVal,:dims,dimPtr]))) + @test length(getproperty(C_dims,dimPtr)) == length(getproperty(dA,dimPtr)) + 1 + @test eltype(getproperty(C_dims,dimPtr)) == eltype(getproperty(dA,dimPtr)) + end + @test getproperty(C_dims,:dims) == (n,n+1) + @testset "similar(A::$h{$elty},Tv::$newElty) $propertyname dims conserved" for propertyname in conserved_dims + @test getproperty(C_dims,propertyname) == getproperty(dA,propertyname) end - for propertyname in structure_conserved_dims - @test length(getproperty(C_dims,propertyname)) ≈ length(getproperty(dA,propertyname)) - @test eltype(getproperty(C_dims,propertyname)) ≈ eltype(getproperty(dA,propertyname)) + + @testset "similar(A::$h{$elty},Tv::$newElty) $propertyname dims structure conserved" for propertyname in structure_conserved_dims + @test length(getproperty(C_dims,propertyname)) == length(getproperty(dA,propertyname)) + @test eltype(getproperty(C_dims,propertyname)) == eltype(getproperty(dA,propertyname)) end - conserved_eltype = Set([nnz,dims,dimPtr,dimVal]) - structure_conserved_eltype = setdiff(properties,union(conserved_eltype,:nzVal)) - @test eltype(getproperty(C_eltype,:nzVal)) ≈ newElty - @test length(getproperty(C_eltype,:nzVal)) ≈ length(getproperty(dA,:nzVal)) - for propertyname in conserved_eltype - @test getproperty(C_dims,propertyname) ≈ getproperty(dA,propertyname) + conserved_eltype = Set([:nnz,:dims,dimPtr,dimVal]) + structure_conserved_eltype = setdiff(properties,union(conserved_eltype,[:nzVal])) + @test eltype(getproperty(C_eltype,:nzVal)) == newElty + @test length(getproperty(C_eltype,:nzVal)) == length(getproperty(dA,:nzVal)) + @testset "similar(A::$h{$elty},Tv::$newElty) $propertyname elty conserved" for propertyname in conserved_eltype + @test getproperty(C_eltype,propertyname) == getproperty(dA,propertyname) end - for propertyname in structure_conserved_eltype - @test length(getproperty(C_dims,propertyname)) ≈ length(getproperty(dA,propertyname)) - @test eltype(getproperty(C_dims,propertyname)) ≈ eltype(getproperty(dA,propertyname)) + + @testset "similar(A::$h{$elty},Tv::$newElty) $propertyname elty structure conserved" for propertyname in structure_conserved_eltype + @test length(getproperty(C_eltype,propertyname)) == length(getproperty(dA,propertyname)) + @test eltype(getproperty(C_eltype,propertyname)) == eltype(getproperty(dA,propertyname)) end - @test eltype(getproperty(C_full,:nzVal)) ≈ newElty - @test length(getproperty(C_full,:nzVal)) ≈ length(getproperty(dA,:nzVal)) - @test length(getproperty(C_full,:dimPtr)) ≈ length(getproperty(dA,:dimPtr)) - @test getproperty(C_dims,:nnz) ≈ getproperty(dA,:nnz) - @test getproperty(C_full,:dims) ≈ (n,n+1) + @testset "similar(A::$h{$elty},Tv::$newElty) full" begin + @test eltype(getproperty(C_full,:nzVal)) == newElty + @test length(getproperty(C_full,:nzVal)) == length(getproperty(dA,:nzVal)) + @test h==CuSparseMatrixCSR ? length(getproperty(C_full,dimPtr)) == length(getproperty(dA,dimPtr)) : length(getproperty(C_full,dimPtr)) == length(getproperty(dA,dimPtr))+1 + @test getproperty(C_dims,:nnz) == getproperty(dA,:nnz) + @test getproperty(C_full,:dims) == (n,n+1) + end end @testset "similar(A::$f($h{$elty}),$newElty)" for elty in typeSet, @@ -74,51 +85,58 @@ n = 10 C_dims = similar(dA,(n,n+1)) C_eltype = similar(dA,newElty) C_full = similar(dA,newElty,(n,n+1)) - @test typeof(C_simple) ≈ typeof(dA) - @test typeof(C_dims) ≈ typeof(dA) - @test typeof(C_eltype) ≈ typeof(dA) - @test typeof(C_full) ≈ typeof(dA) + @test typeof(C_simple) == typeof(parent(dA)) + @test typeof(C_dims) == typeof(parent(dA)) + @test (typeof(C_eltype) == typeof(parent(dA)) && elty==newElty) || ((typeof(C_eltype) != typeof(parent(dA)) && elty!=newElty)) + @test (typeof(C_full) == typeof(parent(dA)) && elty==newElty) || ((typeof(C_full) != typeof(parent(dA)) && elty!=newElty)) - properties = Set([propertynames(dA)]); + properties = Set(propertynames(parent(dA))); - conserved_simple = Set([nnz]) + conserved_simple = Set([:nnz]) structure_conserved_simple = setdiff(properties,conserved_simple); - for propertyname in conserved_simple - @test getproperty(C_simple,propertyname) ≈ getproperty(dA,propertyname) + @testset "similar(A::$f($h{$elty}),$newElty) $propertyname simple conserved" for propertyname in conserved_simple + @test getproperty(C_simple,propertyname) == getproperty(parent(dA),propertyname) end - for propertyname in structure_conserved_simple - @test length(getproperty(C_simple,propertyname)) ≈ length(getproperty(dA,propertyname)) - @test eltype(getproperty(C_simple,propertyname)) ≈ eltype(getproperty(dA,propertyname)) + @testset "similar(A::$f($h{$elty}),$newElty) $propertyname simple structure conserved" for propertyname in structure_conserved_simple + @test length(getproperty(C_simple,propertyname)) == length(getproperty(parent(dA),propertyname)) + @test eltype(getproperty(C_simple,propertyname)) == eltype(getproperty(parent(dA),propertyname)) end - conserved_dims = Set([nnz]) - structure_conserved_dims = setdiff(properties,union(conserved_dims,dimVal,:dims)) - @test getproperty(C_dims,:dims) ≈ (n,n+1) - for propertyname in conserved_dims - @test getproperty(C_dims,propertyname) ≈ getproperty(dA,propertyname) + conserved_dims = Set([:nnz]) + if h==CuSparseMatrixCSR + structure_conserved_dims = setdiff(properties,union(conserved_dims,Set([dimVal,:dims]))) + else #CSC + structure_conserved_dims = setdiff(properties,union(conserved_dims,Set([dimVal,:dims,dimPtr]))) + @test length(getproperty(C_dims,dimPtr)) == length(getproperty(parent(dA),dimPtr)) + 1 + @test eltype(getproperty(C_dims,dimPtr)) == eltype(getproperty(parent(dA),dimPtr)) + end + @test getproperty(C_dims,:dims) == (n,n+1) + @testset "similar(A::$f($h{$elty}),$newElty) $propertyname dims conserved" for propertyname in conserved_dims + @test getproperty(C_dims,propertyname) == getproperty(parent(dA),propertyname) end - for propertyname in structure_conserved_dims - @test length(getproperty(C_dims,propertyname)) ≈ length(getproperty(dA,propertyname)) - @test eltype(getproperty(C_dims,propertyname)) ≈ eltype(getproperty(dA,propertyname)) + @testset "similar(A::$f($h{$elty}),$newElty) $propertyname dims structure conserved" for propertyname in structure_conserved_dims + @test length(getproperty(C_dims,propertyname)) == length(getproperty(parent(dA),propertyname)) + @test eltype(getproperty(C_dims,propertyname)) == eltype(getproperty(parent(dA),propertyname)) end - conserved_eltype = Set([nnz,dims,dimPtr,dimVal]) - structure_conserved_eltype = setdiff(properties,union(conserved_eltype,:nzVal)) - @test eltype(getproperty(C_eltype,:nzVal)) ≈ newElty - @test length(getproperty(C_eltype,:nzVal)) ≈ length(getproperty(dA,:nzVal)) - for propertyname in conserved_eltype - @test getproperty(C_dims,propertyname) ≈ getproperty(dA,propertyname) + conserved_eltype = Set([:nnz,:dims]) + structure_conserved_eltype = setdiff(properties,union(conserved_eltype,[:nzVal])) + @test eltype(getproperty(C_eltype,:nzVal)) == newElty + @test length(getproperty(C_eltype,:nzVal)) == length(getproperty(parent(dA),:nzVal)) + @testset "similar(A::$f($h{$elty}),$newElty) $propertyname elty conserved" for propertyname in conserved_eltype + @test getproperty(C_eltype,propertyname) == getproperty(parent(dA),propertyname) end - for propertyname in structure_conserved_eltype - @test length(getproperty(C_dims,propertyname)) ≈ length(getproperty(dA,propertyname)) - @test eltype(getproperty(C_dims,propertyname)) ≈ eltype(getproperty(dA,propertyname)) + @testset "similar(A::$f($h{$elty}),$newElty) $propertyname elty structure conserved" for propertyname in structure_conserved_eltype + @test length(getproperty(C_eltype,propertyname)) == length(getproperty(parent(dA),propertyname)) + @test eltype(getproperty(C_eltype,propertyname)) == eltype(getproperty(parent(dA),propertyname)) + end + @testset "similar(A::$f($h{$elty}),$newElty) full" begin + @test eltype(getproperty(C_full,:nzVal)) == newElty + @test length(getproperty(C_full,:nzVal)) == length(getproperty(parent(dA),:nzVal)) + @test h==CuSparseMatrixCSR ? length(getproperty(C_full,dimPtr)) == length(getproperty(parent(dA),dimPtr)) : length(getproperty(C_full,dimPtr)) == (length(getproperty(parent(dA),dimPtr))+1) + @test getproperty(C_dims,:nnz) == getproperty(parent(dA),:nnz) + @test getproperty(C_full,:dims) == (n,n+1) end - - @test eltype(getproperty(C_full,:nzVal)) ≈ newElty - @test length(getproperty(C_full,:nzVal)) ≈ length(getproperty(dA,:nzVal)) - @test length(getproperty(C_full,:dimPtr)) ≈ length(getproperty(dA,:dimPtr)) - @test getproperty(C_dims,:nnz) ≈ getproperty(dA,:nnz) - @test getproperty(C_full,:dims) ≈ (n,n+1) end From 4bf8e00ab8f63eb4919afc743adc83fd187ed0cc Mon Sep 17 00:00:00 2001 From: birkmichael <91692218+birkmichael@users.noreply.github.com> Date: Wed, 6 Oct 2021 22:04:18 +0300 Subject: [PATCH 11/21] splatted 1D tuple dims when creating vector --- lib/cusparse/array.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/cusparse/array.jl b/lib/cusparse/array.jl index 711a178fad..473ba18e52 100644 --- a/lib/cusparse/array.jl +++ b/lib/cusparse/array.jl @@ -219,9 +219,9 @@ _cusparsesimilar(S::CuSparseMatrixCSR, ::Type{TvNew}, ::Type{TiNew}, dims::Dims{ CuSparseMatrixCSR(CUDA.fill(one(TiNew), first(dims)+1), similar(colvals(S), TiNew), similar(nonzeros(S), TvNew),dims) # parent method for similar that allocates an empty sparse vector (when new dims are single) _cusparsesimilar(S::CuSparseMatrixCSC, ::Type{TvNew}, ::Type{TiNew}, dims::Dims{1}) where {TvNew,TiNew} = - CuSparseVector(similar(rowvals(S), TiNew, 0), similar(nonzeros(S), TvNew, 0),dims) + CuSparseVector(similar(rowvals(S), TiNew, 0), similar(nonzeros(S), TvNew, 0),dims...) _cusparsesimilar(S::CuSparseMatrixCSR, ::Type{TvNew}, ::Type{TiNew}, dims::Dims{1}) where {TvNew,TiNew} = - CuSparseVector(similar(colvals(S), TiNew, 0), similar(nonzeros(S), TvNew, 0),dims) + CuSparseVector(similar(colvals(S), TiNew, 0), similar(nonzeros(S), TvNew, 0),dims...) """ Utility union type of [`CuSparseMatrixCSC`](@ref), [`CuSparseMatrixCSR`](@ref), for which similar is implemented From 1a097dbc4904179a77d6711815acb6d5e5d0be89 Mon Sep 17 00:00:00 2001 From: birkmichael <91692218+birkmichael@users.noreply.github.com> Date: Sat, 9 Oct 2021 12:54:34 +0300 Subject: [PATCH 12/21] Added method getindex(CuSparseMat,Colon), changed show for vectors --- lib/cusparse/array.jl | 24 ++++++++++++++++++++++-- 1 file changed, 22 insertions(+), 2 deletions(-) diff --git a/lib/cusparse/array.jl b/lib/cusparse/array.jl index 473ba18e52..0061d6bcf6 100644 --- a/lib/cusparse/array.jl +++ b/lib/cusparse/array.jl @@ -316,6 +316,7 @@ Hermitian{T}(Mat::CuSparseMatrix{T}) where T = Hermitian{T,typeof(Mat)}(Mat,'U') # translations Base.getindex(A::AbstractCuSparseVector, ::Colon) = copy(A) +Base.getindex(A::AbstractCuSparseMatrix, ::Colon) = CuSparseVector(A) Base.getindex(A::AbstractCuSparseMatrix, ::Colon, ::Colon) = copy(A) Base.getindex(A::AbstractCuSparseMatrix, i, ::Colon) = getindex(A, i, 1:size(A, 2)) Base.getindex(A::AbstractCuSparseMatrix, ::Colon, i) = getindex(A, 1:size(A, 1), i) @@ -502,8 +503,7 @@ Base.copy(Mat::CuSparseMatrixCOO) = copyto!(similar(Mat), Mat) # input/output -for (gpu, cpu) in [CuSparseVector => SparseVector, - CuSparseMatrixCSC => SparseMatrixCSC, +for (gpu, cpu) in [CuSparseMatrixCSC => SparseMatrixCSC, CuSparseMatrixCSR => SparseMatrixCSC, CuSparseMatrixBSR => SparseMatrixCSC, CuSparseMatrixCOO => SparseMatrixCSC] @@ -529,6 +529,26 @@ for (gpu, cpu) in [CuSparseVector => SparseVector, end +@eval Base.show(io::IOContext, x::CuSparseVector) = + show(io, SparseVector(x)) + +@eval function Base.show(io::IO, mime::MIME"text/plain", S::CuSparseVector) + xnnz = nnz(S) + n = length(S) + print(io, n, "-element ", typeof(S), " with ", xnnz, " stored ", + xnnz == 1 ? "entry" : "entries") + if !(n == 0) + println(io, ":") + io = IOContext(io, :typeinfo => eltype(S)) + if ndims(S) == 1 + show(io, SparseVector(S)) + else + # so that we get the nice Braille pattern + Base.print_array(io, SparseVector(S)) + end + end +end + # interop with device arrays function Adapt.adapt_structure(to::CUDA.Adaptor, x::CuSparseVector) From 38cf586d89f708fde9900ca9c0e1c02783f6fffe Mon Sep 17 00:00:00 2001 From: birkmichael <91692218+birkmichael@users.noreply.github.com> Date: Sat, 9 Oct 2021 12:55:06 +0300 Subject: [PATCH 13/21] Added conversion from matrix to vector --- lib/cusparse/conversions.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/lib/cusparse/conversions.jl b/lib/cusparse/conversions.jl index c016f9e83c..bc9279ae55 100644 --- a/lib/cusparse/conversions.jl +++ b/lib/cusparse/conversions.jl @@ -41,6 +41,8 @@ function SparseArrays.sparse(I::CuVector{Cint}, J::CuVector{Cint}, V::CuVector{T end end +## Matrix to vector, no direct conversion +CuSparseVector(Mat::AbstractCuSparseMatrix) = CuSparseVector(sparse(reshape(Mat,:))) ## CSR to CSC From 575b576cde185c4322de19ddfa89ad9e4dc40e13 Mon Sep 17 00:00:00 2001 From: birkmichael <91692218+birkmichael@users.noreply.github.com> Date: Fri, 22 Oct 2021 02:27:09 +0300 Subject: [PATCH 14/21] Full indexing of CuSparseArrays --- lib/cusparse/array.jl | 68 +++++++++++++++++++++++++++++++++---------- 1 file changed, 52 insertions(+), 16 deletions(-) diff --git a/lib/cusparse/array.jl b/lib/cusparse/array.jl index 8ed05016bb..3f043527d6 100644 --- a/lib/cusparse/array.jl +++ b/lib/cusparse/array.jl @@ -45,6 +45,8 @@ mutable struct CuSparseMatrixCSC{Tv, Ti} <: AbstractCuSparseMatrix{Tv, Ti} end end +CuSparseMatrixCSC(A::CuSparseMatrixCSC) = A + function CUDA.unsafe_free!(xs::CuSparseMatrixCSC) unsafe_free!(xs.colPtr) unsafe_free!(rowvals(xs)) @@ -78,6 +80,8 @@ mutable struct CuSparseMatrixCSR{Tv, Ti} <: AbstractCuSparseMatrix{Tv, Ti} end end +CuSparseMatrixCSR(A::CuSparseMatrixCSR) = A + function CUDA.unsafe_free!(xs::CuSparseMatrixCSR) unsafe_free!(xs.rowPtr) unsafe_free!(xs.colVal) @@ -106,6 +110,8 @@ mutable struct CuSparseMatrixBSR{Tv, Ti} <: AbstractCuSparseMatrix{Tv, Ti} end end +CuSparseMatrixBSR(A::CuSparseMatrixBSR) = A + function CUDA.unsafe_free!(xs::CuSparseMatrixBSR) unsafe_free!(xs.rowPtr) unsafe_free!(xs.colVal) @@ -132,6 +138,8 @@ mutable struct CuSparseMatrixCOO{Tv, Ti} <: AbstractCuSparseMatrix{Tv, Ti} end end +CuSparseMatrixCOO(A::CuSparseMatrixCOO) = A + """ Utility union type of [`CuSparseMatrixCSC`](@ref), [`CuSparseMatrixCSR`](@ref), [`CuSparseMatrixBSR`](@ref), [`CuSparseMatrixCOO`](@ref). @@ -314,24 +322,54 @@ Base.getindex(A::AbstractCuSparseMatrix, i, ::Colon) = getindex(A, i, 1:si Base.getindex(A::AbstractCuSparseMatrix, ::Colon, i) = getindex(A, 1:size(A, 1), i) Base.getindex(A::AbstractCuSparseMatrix, I::Tuple{Integer,Integer}) = getindex(A, I[1], I[2]) -# column slices -function Base.getindex(x::CuSparseMatrixCSC, ::Colon, j::Integer) - checkbounds(x, :, j) - r1 = convert(Int, x.colPtr[j]) - r2 = convert(Int, x.colPtr[j+1]) - 1 - CuSparseVector(rowvals(x)[r1:r2], nonzeros(x)[r1:r2], size(x, 1)) +Base.getindex(A::CuSparseMatrixCSC,I::AbstractVector,J::AbstractVector) = CuSparseMatrixCSC(getindex(CuSparseMatrixCSR(getindex(A,:,J)),I,:)) +Base.getindex(A::CuSparseMatrixCSR,I::AbstractVector,J::AbstractVector) = CuSparseMatrixCSR(getindex(CuSparseMatrixCSC(getindex(A,I,:)),:,J)) + +function Base.getindex(A::CuSparseMatrixCSR,I::AbstractVector,::Colon) + m,n = size(A) + I_sorted = sort(I); + if !(I_sorted[end] <= m && 1 <= I_sorted[1]) + throw(BoundsError()) + end + + rowptr = getrowptr(A) + vals = nonzeros(A) + cols = colvals(A) + + diffVec = vcat(1,diff(rowptr)[I_sorted]) + newrowptr = cumsum(diffVec) + + indices = map((x,y)->range(x,stop=y),rowptr[I_sorted],rowptr[I_sorted.+1].-1) + @CUDA.allowscalar indices = vcat(indices...) #TODO: figure out how to avoid this step + newval = vals[indices] + newcol = cols[indices] + return CuSparseMatrixCSR(newrowptr,newcol,newval,(length(I_sorted),n)) end -function Base.getindex(x::CuSparseMatrixCSR, i::Integer, ::Colon) - checkbounds(x, i, :) - c1 = convert(Int, x.rowPtr[i]) - c2 = convert(Int, x.rowPtr[i+1]) - 1 - CuSparseVector(x.colVal[c1:c2], nonzeros(x)[c1:c2], size(x, 2)) +function Base.getindex(A::CuSparseMatrixCSC,::Colon,I::AbstractVector) + m,n = size(A) + I_sorted = sort(I); + if !(I_sorted[end] <= n && 1 <= I_sorted[1]) + throw(BoundsError()) + end + + colptr = getcolptr(A) + vals = nonzeros(A) + rows = rowvals(A) + + diffVec = vcat(1,diff(colptr)[I_sorted]) + newcolptr = cumsum(diffVec) + + indices = map((x,y)->range(x,stop=y),colptr[I_sorted],colptr[I_sorted.+1].-1) + @CUDA.allowscalar indices = vcat(indices...) #TODO: figure out how to avoid this step + newval = vals[indices] + newrow = rows[indices] + CuSparseMatrixCSC(newcolptr,newrow,newval,(m,length(I))) + return CuSparseMatrixCSC(newcolptr,newrow,newval,(m,length(I))) end -# row slices -Base.getindex(A::CuSparseMatrixCSC, i::Integer, ::Colon) = CuSparseVector(sparse(A[i, 1:end])) # TODO: optimize -Base.getindex(A::CuSparseMatrixCSR, ::Colon, j::Integer) = CuSparseVector(sparse(A[1:end, j])) # TODO: optimize +Base.getindex(x::CuSparseMatrixCSCR, ::Colon, j::Integer) = CuSparseVector(getindex(x,:,[j])) +Base.getindex(x::CuSparseMatrixCSCR, i::Integer ,::Colon) = CuSparseVector(getindex(x,[i],:)) function Base.getindex(A::CuSparseMatrixCSC{T}, i0::Integer, i1::Integer) where T m, n = size(A) @@ -495,7 +533,6 @@ Base.copy(Mat::CuSparseMatrixCOO) = copyto!(similar(Mat), Mat) # input/output - for (gpu, cpu) in [CuSparseVector => SparseVector] @eval function Base.show(io::IO, ::MIME"text/plain", x::$gpu) xnnz = length(nonzeros(x)) @@ -533,7 +570,6 @@ for (gpu, cpu) in [CuSparseMatrixCSC => SparseMatrixCSC, end end - # interop with device arrays function Adapt.adapt_structure(to::CUDA.Adaptor, x::CuSparseVector) From bd16d767be322073037c6d46742648d251474f26 Mon Sep 17 00:00:00 2001 From: birkmichael <91692218+birkmichael@users.noreply.github.com> Date: Fri, 22 Oct 2021 02:28:01 +0300 Subject: [PATCH 15/21] Vectorising matrices --- lib/cusparse/conversions.jl | 43 +++++++++++++++++++++++++------------ 1 file changed, 29 insertions(+), 14 deletions(-) diff --git a/lib/cusparse/conversions.jl b/lib/cusparse/conversions.jl index a27f651119..7b3913060b 100644 --- a/lib/cusparse/conversions.jl +++ b/lib/cusparse/conversions.jl @@ -41,8 +41,23 @@ function SparseArrays.sparse(I::CuVector{Cint}, J::CuVector{Cint}, V::CuVector{T end end -## Matrix to vector, no direct conversion -CuSparseVector(Mat::AbstractCuSparseMatrix) = CuSparseVector(sparse(reshape(Mat,:))) +# CSC to vector +function CuSparseVector(Mat::CuSparseMatrixCSC) + n,m = size(Mat) + dim = n*m + I = 1:m + nzval = nonzeros(Mat) + nzind = rowvals(Mat) + colptr = Array(getcolptr(Mat)) #TODO: figure out how to avoid this step + indices = map((x,y)->range(x,stop=y),colptr[I],colptr[I.+1].-1) + indices = cu.(collect.(indices)) + nzind_array = map(.+,map(inds->nzind[inds],indices) , n.*(0:m-1)) #Vector of cuArrays with indices that should be vectorised into nzind + nzind = vcat(nzind_array...) + return CuSparseVector(nzind,nzval,dim) +end + +# CSR, COO, BSR to vector +CuSparseVector(Mat::AbstractCuSparseMatrix) = CuSparseVector(CuSparseMatrixCSC(Mat)) ## CSR to CSC @@ -221,8 +236,8 @@ for (fname,elty) in ((:cusparseScsr2bsr, :Float32), mb = div((m + blockDim - 1),blockDim) nb = div((n + blockDim - 1),blockDim) bsrRowPtr = CUDA.zeros(Cint,mb + 1) - cudesca = CuMatrixDescriptor('G', 'L', 'N', inda) - cudescc = CuMatrixDescriptor('G', 'L', 'N', indc) + cudesca = CuMatrixDescriptor(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, inda) + cudescc = CuMatrixDescriptor(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, indc) cusparseXcsr2bsrNnz(handle(), dir, m, n, cudesca, csr.rowPtr, csr.colVal, blockDim, cudescc, bsrRowPtr, nnz_ref) bsrNzVal = CUDA.zeros($elty, nnz_ref[] * blockDim * blockDim ) @@ -247,8 +262,8 @@ for (fname,elty) in ((:cusparseSbsr2csr, :Float32), mb = div(m,bsr.blockDim) nb = div(n,bsr.blockDim) nnzVal = nnz(bsr) * bsr.blockDim * bsr.blockDim - cudesca = CuMatrixDescriptor('G', 'L', 'N', inda) - cudescc = CuMatrixDescriptor('G', 'L', 'N', indc) + cudesca = CuMatrixDescriptor(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, inda) + cudescc = CuMatrixDescriptor(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, indc) csrRowPtr = CUDA.zeros(Cint, m + 1) csrColInd = CUDA.zeros(Cint, nnzVal) csrNzVal = CUDA.zeros($elty, nnzVal) @@ -313,7 +328,7 @@ for (cname,rname,elty) in ((:cusparseScsc2dense, :cusparseScsr2dense, :Float32), end return denseA else - cudesc = CuMatrixDescriptor('G', 'L', 'N', ind) + cudesc = CuMatrixDescriptor(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, ind) lda = max(1,stride(denseA,2)) $rname(handle(), m, n, cudesc, nonzeros(csr), csr.rowPtr, csr.colVal, denseA, lda) @@ -340,7 +355,7 @@ for (cname,rname,elty) in ((:cusparseScsc2dense, :cusparseScsr2dense, :Float32), return denseA else lda = max(1,stride(denseA,2)) - cudesc = CuMatrixDescriptor('G', 'L', 'N', ind) + cudesc = CuMatrixDescriptor(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, ind) $cname(handle(), m, n, cudesc, nonzeros(csc), rowvals(csc), csc.colPtr, denseA, lda) return denseA @@ -415,9 +430,9 @@ for (nname,cname,rname,elty) in ((:cusparseSnnz, :cusparseSdense2csc, :cusparseS function CuSparseMatrixCSR(A::CuMatrix{$elty}; ind::SparseChar='O') m,n = size(A) lda = max(1, stride(A,2)) - cudesc = CuMatrixDescriptor('G', - 'L', - 'N', ind) + cudesc = CuMatrixDescriptor(CUSPARSE_MATRIX_TYPE_GENERAL, + CUSPARSE_FILL_MODE_LOWER, + CUSPARSE_DIAG_TYPE_NON_UNIT, ind) nnzRowCol = CUDA.zeros(Cint, m) nnzTotal = Ref{Cint}(1) $nname(handle(), @@ -435,9 +450,9 @@ for (nname,cname,rname,elty) in ((:cusparseSnnz, :cusparseSdense2csc, :cusparseS function CuSparseMatrixCSC(A::CuMatrix{$elty}; ind::SparseChar='O') m,n = size(A) lda = max(1, stride(A,2)) - cudesc = CuMatrixDescriptor('G', - 'L', - 'N', ind) + cudesc = CuMatrixDescriptor(CUSPARSE_MATRIX_TYPE_GENERAL, + CUSPARSE_FILL_MODE_LOWER, + CUSPARSE_DIAG_TYPE_NON_UNIT, ind) nnzRowCol = CUDA.zeros(Cint, n) nnzTotal = Ref{Cint}(1) $nname(handle(), From 69e69d6dd07dc991a7ea0be114d4a9f4a3feb37f Mon Sep 17 00:00:00 2001 From: birkmichael <91692218+birkmichael@users.noreply.github.com> Date: Fri, 22 Oct 2021 13:37:02 +0300 Subject: [PATCH 16/21] Added vectorising of CuSparseMatrices --- lib/cusparse/conversions.jl | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/lib/cusparse/conversions.jl b/lib/cusparse/conversions.jl index 7b3913060b..21fa970dbe 100644 --- a/lib/cusparse/conversions.jl +++ b/lib/cusparse/conversions.jl @@ -236,8 +236,8 @@ for (fname,elty) in ((:cusparseScsr2bsr, :Float32), mb = div((m + blockDim - 1),blockDim) nb = div((n + blockDim - 1),blockDim) bsrRowPtr = CUDA.zeros(Cint,mb + 1) - cudesca = CuMatrixDescriptor(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, inda) - cudescc = CuMatrixDescriptor(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, indc) + cudesca = CuMatrixDescriptor('G', 'L', 'N', inda) + cudescc = CuMatrixDescriptor('G', 'L', 'N', indc) cusparseXcsr2bsrNnz(handle(), dir, m, n, cudesca, csr.rowPtr, csr.colVal, blockDim, cudescc, bsrRowPtr, nnz_ref) bsrNzVal = CUDA.zeros($elty, nnz_ref[] * blockDim * blockDim ) @@ -262,8 +262,8 @@ for (fname,elty) in ((:cusparseSbsr2csr, :Float32), mb = div(m,bsr.blockDim) nb = div(n,bsr.blockDim) nnzVal = nnz(bsr) * bsr.blockDim * bsr.blockDim - cudesca = CuMatrixDescriptor(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, inda) - cudescc = CuMatrixDescriptor(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, indc) + cudesca = CuMatrixDescriptor('G', 'L', 'N', inda) + cudescc = CuMatrixDescriptor('G', 'L', 'N', indc) csrRowPtr = CUDA.zeros(Cint, m + 1) csrColInd = CUDA.zeros(Cint, nnzVal) csrNzVal = CUDA.zeros($elty, nnzVal) @@ -313,7 +313,7 @@ for (cname,rname,elty) in ((:cusparseScsc2dense, :cusparseScsr2dense, :Float32), m,n = csr.dims denseA = CUDA.zeros($elty,m,n) if version() >= v"11.3.0" # CUSPARSE version from CUDA release notes - desc_csr = CuSparseMatrixDescriptor(csr) + desc_csr = CuSparseMatrixDescriptor(csr, ind) desc_dense = CuDenseMatrixDescriptor(denseA) function bufferSize() @@ -328,7 +328,7 @@ for (cname,rname,elty) in ((:cusparseScsc2dense, :cusparseScsr2dense, :Float32), end return denseA else - cudesc = CuMatrixDescriptor(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, ind) + cudesc = CuMatrixDescriptor('G', 'L', 'N', ind) lda = max(1,stride(denseA,2)) $rname(handle(), m, n, cudesc, nonzeros(csr), csr.rowPtr, csr.colVal, denseA, lda) @@ -339,7 +339,7 @@ for (cname,rname,elty) in ((:cusparseScsc2dense, :cusparseScsr2dense, :Float32), m,n = csc.dims denseA = CUDA.zeros($elty,m,n) if version() >= v"11.3.0" # CUSPARSE version from CUDA release notes - desc_csc = CuSparseMatrixDescriptor(csc; convert=false) + desc_csc = CuSparseMatrixDescriptor(csc, ind; convert=false) desc_dense = CuDenseMatrixDescriptor(denseA) function bufferSize() @@ -355,7 +355,7 @@ for (cname,rname,elty) in ((:cusparseScsc2dense, :cusparseScsr2dense, :Float32), return denseA else lda = max(1,stride(denseA,2)) - cudesc = CuMatrixDescriptor(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, ind) + cudesc = CuMatrixDescriptor('G', 'L', 'N', ind) $cname(handle(), m, n, cudesc, nonzeros(csc), rowvals(csc), csc.colPtr, denseA, lda) return denseA @@ -371,7 +371,7 @@ for (elty, welty) in ((:Float16, :Float32), m,n = csr.dims denseA = CUDA.zeros($elty,m,n) if version() >= v"11.3.0" # CUSPARSE version from CUDA release notes - desc_csr = CuSparseMatrixDescriptor(csr) + desc_csr = CuSparseMatrixDescriptor(csr, ind) desc_dense = CuDenseMatrixDescriptor(denseA) function bufferSize() @@ -396,7 +396,7 @@ for (elty, welty) in ((:Float16, :Float32), m,n = csc.dims denseA = CUDA.zeros($elty,m,n) if version() >= v"11.3.0" # CUSPARSE version from CUDA release notes - desc_csc = CuSparseMatrixDescriptor(csc; convert=false) + desc_csc = CuSparseMatrixDescriptor(csc, ind; convert=false) desc_dense = CuDenseMatrixDescriptor(denseA) function bufferSize() @@ -430,9 +430,9 @@ for (nname,cname,rname,elty) in ((:cusparseSnnz, :cusparseSdense2csc, :cusparseS function CuSparseMatrixCSR(A::CuMatrix{$elty}; ind::SparseChar='O') m,n = size(A) lda = max(1, stride(A,2)) - cudesc = CuMatrixDescriptor(CUSPARSE_MATRIX_TYPE_GENERAL, - CUSPARSE_FILL_MODE_LOWER, - CUSPARSE_DIAG_TYPE_NON_UNIT, ind) + cudesc = CuMatrixDescriptor('G', + 'L', + 'N', ind) nnzRowCol = CUDA.zeros(Cint, m) nnzTotal = Ref{Cint}(1) $nname(handle(), @@ -450,9 +450,9 @@ for (nname,cname,rname,elty) in ((:cusparseSnnz, :cusparseSdense2csc, :cusparseS function CuSparseMatrixCSC(A::CuMatrix{$elty}; ind::SparseChar='O') m,n = size(A) lda = max(1, stride(A,2)) - cudesc = CuMatrixDescriptor(CUSPARSE_MATRIX_TYPE_GENERAL, - CUSPARSE_FILL_MODE_LOWER, - CUSPARSE_DIAG_TYPE_NON_UNIT, ind) + cudesc = CuMatrixDescriptor('G', + 'L', + 'N', ind) nnzRowCol = CUDA.zeros(Cint, n) nnzTotal = Ref{Cint}(1) $nname(handle(), From 9a7196f622510dc29ed8ee654490685e56ad0135 Mon Sep 17 00:00:00 2001 From: birkmichael <91692218+birkmichael@users.noreply.github.com> Date: Fri, 22 Oct 2021 16:06:18 +0300 Subject: [PATCH 17/21] Update array.jl --- lib/cusparse/array.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/lib/cusparse/array.jl b/lib/cusparse/array.jl index 3f043527d6..773f6060fd 100644 --- a/lib/cusparse/array.jl +++ b/lib/cusparse/array.jl @@ -370,6 +370,8 @@ end Base.getindex(x::CuSparseMatrixCSCR, ::Colon, j::Integer) = CuSparseVector(getindex(x,:,[j])) Base.getindex(x::CuSparseMatrixCSCR, i::Integer ,::Colon) = CuSparseVector(getindex(x,[i],:)) +Base.getindex(x::CuSparseMatrixCSCR, I::AbstractVector, j::Integer) = CuSparseVector(getindex(x, I, [j])) +Base.getindex(x::CuSparseMatrixCSCR, i::Integer, J::AbstractVector) = CuSparseVector(getindex(x, [i], J)) function Base.getindex(A::CuSparseMatrixCSC{T}, i0::Integer, i1::Integer) where T m, n = size(A) From f90a4d66564e1978a8241a1fac36921f90a8ad4f Mon Sep 17 00:00:00 2001 From: birkmichael <91692218+birkmichael@users.noreply.github.com> Date: Fri, 22 Oct 2021 16:12:29 +0300 Subject: [PATCH 18/21] Update conversions.jl --- lib/cusparse/conversions.jl | 25 +++++++++++++++++++++---- 1 file changed, 21 insertions(+), 4 deletions(-) diff --git a/lib/cusparse/conversions.jl b/lib/cusparse/conversions.jl index 81378de33c..915b9e68c1 100644 --- a/lib/cusparse/conversions.jl +++ b/lib/cusparse/conversions.jl @@ -41,6 +41,23 @@ function SparseArrays.sparse(I::CuVector{Cint}, J::CuVector{Cint}, V::CuVector{T end end +# CSC to vector +function CuSparseVector(Mat::CuSparseMatrixCSC) + n,m = size(Mat) + dim = n*m + I = 1:m + nzval = nonzeros(Mat) + nzind = rowvals(Mat) + colptr = Array(getcolptr(Mat)) #TODO: figure out how to avoid this step + indices = map((x,y)->range(x,stop=y),colptr[I],colptr[I.+1].-1) + indices = cu.(collect.(indices)) + nzind_array = map(.+,map(inds->nzind[inds],indices) , n.*(0:m-1)) #Vector of cuArrays with indices that should be vectorised into nzind + nzind = vcat(nzind_array...) + return CuSparseVector(nzind,nzval,dim) +end + +# CSR, COO, BSR to vector +CuSparseVector(Mat::AbstractCuSparseMatrix) = CuSparseVector(CuSparseMatrixCSC(Mat)) ## CSR to CSC @@ -296,7 +313,7 @@ for (cname,rname,elty) in ((:cusparseScsc2dense, :cusparseScsr2dense, :Float32), m,n = csr.dims denseA = CUDA.zeros($elty,m,n) if version() >= v"11.3.0" # CUSPARSE version from CUDA release notes - desc_csr = CuSparseMatrixDescriptor(csr, ind) + desc_csr = CuSparseMatrixDescriptor(csr) desc_dense = CuDenseMatrixDescriptor(denseA) function bufferSize() @@ -322,7 +339,7 @@ for (cname,rname,elty) in ((:cusparseScsc2dense, :cusparseScsr2dense, :Float32), m,n = csc.dims denseA = CUDA.zeros($elty,m,n) if version() >= v"11.3.0" # CUSPARSE version from CUDA release notes - desc_csc = CuSparseMatrixDescriptor(csc, ind; convert=false) + desc_csc = CuSparseMatrixDescriptor(csc; convert=false) desc_dense = CuDenseMatrixDescriptor(denseA) function bufferSize() @@ -354,7 +371,7 @@ for (elty, welty) in ((:Float16, :Float32), m,n = csr.dims denseA = CUDA.zeros($elty,m,n) if version() >= v"11.3.0" # CUSPARSE version from CUDA release notes - desc_csr = CuSparseMatrixDescriptor(csr, ind) + desc_csr = CuSparseMatrixDescriptor(csr) desc_dense = CuDenseMatrixDescriptor(denseA) function bufferSize() @@ -379,7 +396,7 @@ for (elty, welty) in ((:Float16, :Float32), m,n = csc.dims denseA = CUDA.zeros($elty,m,n) if version() >= v"11.3.0" # CUSPARSE version from CUDA release notes - desc_csc = CuSparseMatrixDescriptor(csc, ind; convert=false) + desc_csc = CuSparseMatrixDescriptor(csc; convert=false) desc_dense = CuDenseMatrixDescriptor(denseA) function bufferSize() From 0e45e485877081476073b1098db70cd7c553df42 Mon Sep 17 00:00:00 2001 From: birkmichael <91692218+birkmichael@users.noreply.github.com> Date: Fri, 22 Oct 2021 16:14:27 +0300 Subject: [PATCH 19/21] Update conversions.jl --- lib/cusparse/conversions.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/lib/cusparse/conversions.jl b/lib/cusparse/conversions.jl index 21fa970dbe..915b9e68c1 100644 --- a/lib/cusparse/conversions.jl +++ b/lib/cusparse/conversions.jl @@ -313,7 +313,7 @@ for (cname,rname,elty) in ((:cusparseScsc2dense, :cusparseScsr2dense, :Float32), m,n = csr.dims denseA = CUDA.zeros($elty,m,n) if version() >= v"11.3.0" # CUSPARSE version from CUDA release notes - desc_csr = CuSparseMatrixDescriptor(csr, ind) + desc_csr = CuSparseMatrixDescriptor(csr) desc_dense = CuDenseMatrixDescriptor(denseA) function bufferSize() @@ -339,7 +339,7 @@ for (cname,rname,elty) in ((:cusparseScsc2dense, :cusparseScsr2dense, :Float32), m,n = csc.dims denseA = CUDA.zeros($elty,m,n) if version() >= v"11.3.0" # CUSPARSE version from CUDA release notes - desc_csc = CuSparseMatrixDescriptor(csc, ind; convert=false) + desc_csc = CuSparseMatrixDescriptor(csc; convert=false) desc_dense = CuDenseMatrixDescriptor(denseA) function bufferSize() @@ -371,7 +371,7 @@ for (elty, welty) in ((:Float16, :Float32), m,n = csr.dims denseA = CUDA.zeros($elty,m,n) if version() >= v"11.3.0" # CUSPARSE version from CUDA release notes - desc_csr = CuSparseMatrixDescriptor(csr, ind) + desc_csr = CuSparseMatrixDescriptor(csr) desc_dense = CuDenseMatrixDescriptor(denseA) function bufferSize() @@ -396,7 +396,7 @@ for (elty, welty) in ((:Float16, :Float32), m,n = csc.dims denseA = CUDA.zeros($elty,m,n) if version() >= v"11.3.0" # CUSPARSE version from CUDA release notes - desc_csc = CuSparseMatrixDescriptor(csc, ind; convert=false) + desc_csc = CuSparseMatrixDescriptor(csc; convert=false) desc_dense = CuDenseMatrixDescriptor(denseA) function bufferSize() From e62cbb6ef3c47f97c405a886a7104d07614fdac0 Mon Sep 17 00:00:00 2001 From: birkmichael <91692218+birkmichael@users.noreply.github.com> Date: Fri, 22 Oct 2021 16:19:38 +0300 Subject: [PATCH 20/21] changed to match updates in master --- lib/cusparse/conversions.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/lib/cusparse/conversions.jl b/lib/cusparse/conversions.jl index 915b9e68c1..61a65bf56c 100644 --- a/lib/cusparse/conversions.jl +++ b/lib/cusparse/conversions.jl @@ -313,7 +313,7 @@ for (cname,rname,elty) in ((:cusparseScsc2dense, :cusparseScsr2dense, :Float32), m,n = csr.dims denseA = CUDA.zeros($elty,m,n) if version() >= v"11.3.0" # CUSPARSE version from CUDA release notes - desc_csr = CuSparseMatrixDescriptor(csr) + desc_csr = CuSparseMatrixDescriptor(csr ,ind) desc_dense = CuDenseMatrixDescriptor(denseA) function bufferSize() @@ -339,7 +339,7 @@ for (cname,rname,elty) in ((:cusparseScsc2dense, :cusparseScsr2dense, :Float32), m,n = csc.dims denseA = CUDA.zeros($elty,m,n) if version() >= v"11.3.0" # CUSPARSE version from CUDA release notes - desc_csc = CuSparseMatrixDescriptor(csc; convert=false) + desc_csc = CuSparseMatrixDescriptor(csc, ind; convert=false) desc_dense = CuDenseMatrixDescriptor(denseA) function bufferSize() @@ -371,7 +371,7 @@ for (elty, welty) in ((:Float16, :Float32), m,n = csr.dims denseA = CUDA.zeros($elty,m,n) if version() >= v"11.3.0" # CUSPARSE version from CUDA release notes - desc_csr = CuSparseMatrixDescriptor(csr) + desc_csr = CuSparseMatrixDescriptor(csr, ind) desc_dense = CuDenseMatrixDescriptor(denseA) function bufferSize() @@ -396,7 +396,7 @@ for (elty, welty) in ((:Float16, :Float32), m,n = csc.dims denseA = CUDA.zeros($elty,m,n) if version() >= v"11.3.0" # CUSPARSE version from CUDA release notes - desc_csc = CuSparseMatrixDescriptor(csc; convert=false) + desc_csc = CuSparseMatrixDescriptor(csc, ind; convert=false) desc_dense = CuDenseMatrixDescriptor(denseA) function bufferSize() From 509d10a907555a0f439194971267d0c5d43825d4 Mon Sep 17 00:00:00 2001 From: birkmichael <91692218+birkmichael@users.noreply.github.com> Date: Fri, 22 Oct 2021 16:20:28 +0300 Subject: [PATCH 21/21] Update conversions.jl --- lib/cusparse/conversions.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/cusparse/conversions.jl b/lib/cusparse/conversions.jl index 61a65bf56c..21fa970dbe 100644 --- a/lib/cusparse/conversions.jl +++ b/lib/cusparse/conversions.jl @@ -313,7 +313,7 @@ for (cname,rname,elty) in ((:cusparseScsc2dense, :cusparseScsr2dense, :Float32), m,n = csr.dims denseA = CUDA.zeros($elty,m,n) if version() >= v"11.3.0" # CUSPARSE version from CUDA release notes - desc_csr = CuSparseMatrixDescriptor(csr ,ind) + desc_csr = CuSparseMatrixDescriptor(csr, ind) desc_dense = CuDenseMatrixDescriptor(denseA) function bufferSize()