# Looking for code of sparse matrix \* dense matrix
**URL:** https://discourse.julialang.org/t/looking-for-code-of-sparse-matrix-dense-matrix/60283
**Category:** Internals & Design
**Created:** [April 30, 2021, 1:18am UTC](https://discourse.julialang.org/t/looking-for-code-of-sparse-matrix-dense-matrix/60283 "2021-04-30T01:18:31Z")
**Posts on this page:** 8
**Page:** 1
### Author:  [@Lecrapouille](https://discourse.julialang.org/u/Lecrapouille)
#### Post date: [April 30, 2021, 1:18am UTC](https://discourse.julialang.org/t/looking-for-code-of-sparse-matrix-dense-matrix/60283/1 "2021-04-30T01:18:32Z")
Hi ! I may have found an issue in Julia 1.6 when doing matrix product between a sparse matrix and a full matrix when function zero() and one() and operators + and \* have been changed for a particular algebra. I’m searching inside the code source of Julia in julia/stdlib/SparseArrays/src/sparsematrix.jl for this specific matrix product but I did not find it. Can somebody show me the function name. I’ll analyze the code before opening an issue.
Thanks in advance !
```julia
julia> A
2×2 SparseMatrixCSC{MP{Float64}, Int64} with 4 stored entries:
# Wrong answer (same for full(A)*A)
julia> A*full(A)
# Good answer
julia> A*A
```
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [April 30, 2021, 1:47am UTC](https://discourse.julialang.org/t/looking-for-code-of-sparse-matrix-dense-matrix/60283/2 "2021-04-30T01:47:46Z")
Can you give a MWE (definition of `MP` and associated methods?). You can find the method signature with `@which`, which results in the method `*(A::SparseMatrixCSCUnion{TA}, B::AdjOrTransDenseMatrix) where {TA}` in `stdlib/v1.7/SparseArrays/src/linalg.jl`
---
### Author:  [@ettersi](https://discourse.julialang.org/u/ettersi)
#### Post date: [April 30, 2021, 9:13am UTC](https://discourse.julialang.org/t/looking-for-code-of-sparse-matrix-dense-matrix/60283/3 "2021-04-30T09:13:05Z")
A debugger is also very useful for these purposes. For example, in VS Code type `@enter A*B`, then step inside functions until you find the code that you are interested in.
---
### Author:  [@Lecrapouille](https://discourse.julialang.org/u/Lecrapouille)
#### Post date: [April 30, 2021, 6:57pm UTC](https://discourse.julialang.org/t/looking-for-code-of-sparse-matrix-dense-matrix/60283/4 "2021-04-30T18:57:08Z")
Hi ! Thanks for the link I’ll watch.
Here is the miminal code extracted from my project. I’ll try with the debugger 🙂 but it seems to me this problem has been introduced since Julia 1.4.2 because I tested on versions 0.7 1.2.x 1.3.x and got the correct answers.
MP means Max-Plus: a tropical algebra (see here for more information [http://www.stanczyk.pro/mpa/download.php?id=max-plus-1.7&ex=pdf](http://www.stanczyk.pro/mpa/download.php?id=max-plus-1.7&ex=pdf) )
This code redefines operators + as operator max() and operator \* as operator \*. We define neutral and absorbing elements zeros() as -Inf and one() as 0. I did not include conversion and promotion (not used in this basic example). The magic in Julia is that all complex algebra functions are automatically working … when there is no regression introduced by newer version 🙂
```julia
using LinearAlgebra, SparseArrays
struct MP
v::Float64
end
Base.zero(::Type{MP}) = MP(typemin(Float64))
Base.zero(x::MP) = zero(typeof(x))
Base.one(::Type{MP}) = MP(zero(Float64))
Base.one(x::MP) = one(typeof(x))
Base.:(+)(x::MP, y::MP) = MP(max(x.v, y.v))
Base.:(+)(x::MP, y::Real) = MP(max(x.v, y))
Base.:(+)(x::Real, y::MP) = MP(max(x, y.v))
Base.:(*)(x::MP, y::MP) = MP(x.v + y.v)
Base.:(*)(x::MP, y::Real) = MP(x.v + y)
Base.:(*)(x::Real, y::MP) = MP(x + y.v)
#Not needed for this example:
#MP(x::MP) = MP(x.v)
#Base.promote_rule(::Type{MP}, ::Type{U}) where {U} = MP
#Base.convert(::MP, x::Number) = MP(x)
```
This code gives the good result:
```julia
julia> A=[MP(1) MP(2); MP(3) MP(4)]
julia> A*A
2×2 Array{MP,2}:
MP(5.0) MP(6.0)
MP(7.0) MP(8.0)
```
This one does not give the correct answer since Julia 1.4.2 (1.5.x, 1.6.0):
```julia
julia> sparse(A)*A
2×2 Matrix{MP}:
MP(6.0) MP(7.0)
MP(8.0) MP(9.0)
```
When commenting
```julia
#Base.:(+)(x::MP, y::Real) = MP(max(x.v, y))
#Base.:(+)(x::Real, y::MP) = MP(max(x, y.v))
#Base.:(*)(x::MP, y::Real) = MP(x.v + y)
#Base.:(*)(x::Real, y::MP) = MP(x + y.v)
```
I’ll have this interesting error:
```julia
ERROR: MethodError: no method matching *(::MP, ::Bool)
Closest candidates are:
*(::Any, ::Any, ::Any, ::Any...) at operators.jl:560
*(::StridedArray{P, N} where N, ::Real) where P<:Dates.Period at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/Dates/src/deprecated.jl:44
*(::Union{SparseVector{Tv, Ti}, SubArray{Tv, 1, var"#s832", Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, false} where var"#s832"<:SparseArrays.AbstractSparseMatrixCSC{Tv, Ti}, SubArray{Tv, 1, var"#s832", Tuple{Base.Slice{Base.OneTo{Int64}}}, false} where var"#s832"<:AbstractSparseVector{Tv, Ti}} where {Tv, Ti}, ::Number) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/sparsevector.jl:1448
...
Stacktrace:
[1] mul!(C::Matrix{MP}, A::SparseMatrixCSC{MP, Int64}, B::Matrix{MP}, α::Bool, β::Bool)
@ SparseArrays /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/linalg.jl:41
[2] *(A::SparseMatrixCSC{MP, Int64}, B::Matrix{MP})
@ SparseArrays /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/linalg.jl:53
[3] top-level scope
@ REPL[12]:1
```
This is not good since we are not supposed to call `*(::MP, ::Bool)` bool is not for Max-Plus algebra. Note that we may want to write `MP(1.0) + 3.0` just for the user to avoid typing too many characters 3.0 meaning MP(3.0)) but `MP(1.0) + false` has nonsense.
By the way this is a second regression I found. Identity matrix that used to be created with the `eye()` function was made like this:
```julia
[ one() zero() zero()
zero() one() zero()
zero() zero() one()]
```
Since `eye()` no longer exist, replaced in Julia 1.x (I guess) by `I` which is not working because AGAIN with introduction boolean instead of zero() and one()
```julia
julia> I
UniformScaling{Bool}
true*I
```
and therefore Identity are no longer working in Max-Plus. I succeeded to make a workaround in my package but I have to find a real fix. See [Matrix identity not calling one() ? · Issue #33036 · JuliaLang/julia · GitHub](https://github.com/JuliaLang/julia/issues/33036)
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [April 30, 2021, 7:24pm UTC](https://discourse.julialang.org/t/looking-for-code-of-sparse-matrix-dense-matrix/60283/5 "2021-04-30T19:24:32Z")
Ah, the problem is that sparse matrices assumes that `false` will work as a strong zero (`false+x=x`, and `false*x=false`) for your type. It would be a bit weird, but defining
`MP(x::Bool)= x ? one(MP) : zero(MP)`
should fix it.
---
### Author:  [@Lecrapouille](https://discourse.julialang.org/u/Lecrapouille)
#### Post date: [April 30, 2021, 7:55pm UTC](https://discourse.julialang.org/t/looking-for-code-of-sparse-matrix-dense-matrix/60283/6 "2021-04-30T19:55:29Z")
It does not seem to fix. Tomorrow, I’ll try to diff the julia code between 1.3.x and 1.4.x and see if something odd has changed.
---
### Author:  [@ettersi](https://discourse.julialang.org/u/ettersi)
#### Post date: [May 4, 2021, 5:41am UTC](https://discourse.julialang.org/t/looking-for-code-of-sparse-matrix-dense-matrix/60283/7 "2021-05-04T05:41:34Z")
Regarding your specific question (locate the source code for the dense times sparse matrix product):
[https://github.com/JuliaLang/julia/blob/e3adc399e6bf69d0801c243f254ff8624f50d3c4/stdlib/SparseArrays/src/linalg.jl#L87-L115](https://github.com/JuliaLang/julia/blob/e3adc399e6bf69d0801c243f254ff8624f50d3c4/stdlib/SparseArrays/src/linalg.jl#L87-L115)
Regarding the broader issue (implement the max-plus algebra): I would recommend against implementing `Number <-> MP` conversion and mixed arithmetic. On a more philosophical level, the problem with mixed arithmetic is that it is not clear whether `MP(1) + 2` should be evaluated according to standard or max-plus arithmetic. In more practical terms, the issue is that the Julia standard library implicitly assumes standard arithmetic in at least a couple of places. This leads to silent errors if you do implement mixed arithmetic, whereas if you don’t then you get an error message which is much easier to detect and address.
---
### Author:  [@Lecrapouille](https://discourse.julialang.org/u/Lecrapouille)
#### Post date: [May 5, 2021, 11:03pm UTC](https://discourse.julialang.org/t/looking-for-code-of-sparse-matrix-dense-matrix/60283/8 "2021-05-05T23:03:04Z")
@attersi Thank you for your search, I did not get time to dive into code since hard for me to understand and missing comments.
For the moment since sparse \* sparse is working I made a fallback for dense matrix \* sparse where I convert the dense to sparse. This is allowing me to get the same results as in previous Julia versions. Once my work done i’ll get back investigating on this regression.
Concerning Number \<-\> MP: contamination is fine since working with MP number make you dive in (max,+) algebra and no more in classic algebra. For classic addition: MP.v + scalar is okay