Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Multiplying triangular sparse matrix with sparse matrix fails to return sparse matrix #35610

Closed
KlausC opened this issue Apr 27, 2020 · 2 comments · Fixed by #35659
Closed

Multiplying triangular sparse matrix with sparse matrix fails to return sparse matrix #35610

KlausC opened this issue Apr 27, 2020 · 2 comments · Fixed by #35659
Labels

Comments

@KlausC
Copy link
Contributor

KlausC commented Apr 27, 2020

Multiplying an UpperTriangular{SparseMatrixCSC} with a SparseMatrixCSC should return a sparse matrix. It returns a dense matrix as of version 1.2.
Computational effort should be approximately 1/2 of sparse multiplication of SparseMatrixCSC.

julia> VERSION
v"1.4.1" # same for v"1.2.0"
julia> X = UpperTriangular(sprand(100, 100, 10*1e-2));
julia> Y = sprand(100, 10, 10*1e-2);
julia> typeof(X * Y)
Array{Float64,2}

With version 1, we have the same error as in #35609

julia> VERSION
v"1.0.5"
julia> X = UnitUpperTriangular(sprand(100, 100, 10*1e-2));
julia> Y = sprand(100, 10, 10*1e-2);
julia> X * Y;
ERROR: MethodError: no method matching lmul!(::UnitUpperTriangular{Float64,SparseMatrixCSC{Float64,Int64}}, ::SparseMatrixCSC{Float64,Int64})
@Micket
Copy link
Contributor

Micket commented Apr 30, 2020

(\)(L::TriangularSparse, B::AbstractSparseMatrixCSC) = ldiv!(L, Array(B))
(*)(L::TriangularSparse, B::AbstractSparseMatrixCSC) = lmul!(L, Array(B))

Looks to me that there was some type of attempt to treat * and \ similarly, though, I can't think of any case where you actually would want to do that.

Y*X actually does the right thing there;

*(A::AbstractSparseMatrixCSC{TA}, B::AdjOrTransStridedOrTriangularMatrix{Tx}) where {TA,Tx} =
(T = promote_op(matprod, TA, Tx); mul!(similar(B, T, (size(A, 1), size(B, 2))), A, B, true, false))

and, as far as I can see, this methods exists and does the right thing already:
*(X::AdjOrTransStridedOrTriangularMatrix, A::AbstractSparseMatrixCSC{TvA}) where {TvA} =
(T = promote_op(matprod, eltype(X), TvA); mul!(similar(X, T, (size(X, 1), size(A, 2))), X, A, true, false))

and already does the correct thing.

So, i think this should simply be removed:

(*)(L::TriangularSparse, B::AbstractSparseMatrixCSC) = lmul!(L, Array(B))

(edit: I totally missed that KlausC already had a PR going)

@KlausC
Copy link
Contributor Author

KlausC commented May 1, 2020

Y*X actually does the right thing there;

The point of T * A is the return type, which should be sparse, not dense.

as far as I can see, this methods exists and does the right thing already

The problem with A * T is the execution efficiency. mul! is using a generic implementation here, which does not take much advantage of the sparsity.

The pending PR fixed both cases using the famous Gustafson multiplication algorithm.

StefanKarpinski pushed a commit that referenced this issue Jun 1, 2020
…5659)

* multiplication of sparse triangular matrices
* removed disambiguity and improved test cases
* test cases for `nnz, nzrange, rowvals, nonzeros`
simeonschaub pushed a commit to simeonschaub/julia that referenced this issue Aug 11, 2020
…g#35610 JuliaLang#35642 (JuliaLang#35659)

* multiplication of sparse triangular matrices
* removed disambiguity and improved test cases
* test cases for `nnz, nzrange, rowvals, nonzeros`
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

Successfully merging a pull request may close this issue.

4 participants