Huge v1.8 regression in sparse matrix division with diagonal matrix

[wsshin]

(The original issue was posted at JuliaLang/LinearAlgebra.jl#970, but I was recommended to re-post it here.)

The two sparse matrices A and B are constructed by the finite difference method for partial differential equations. A corresponds to a nearly diagonal matrix representing material parameters; B corresponds to a curl operator in a 2D Cartesian space. Here is how they look like:

julia> A
160000×160000 SparseMatrixCSC{ComplexF64, Int64} with 160000 stored entries:
⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄

julia> B
160000×80000 SparseMatrixCSC{ComplexF64, Int64} with 320000 stored entries:
⢳⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⢳⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⢳⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⢳⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⢳⡀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⢳⡀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠈⢷⡀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⢷⡀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠈⢷⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠈⢧⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢧⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢧⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢧⠀
⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢧

Now, Here is the timing of first execution of A \ B in Julia 1.8. It takes 557 seconds:

julia> VERSION
v"1.8.3-pre.2"

julia> @timev A \ B
557.904421 seconds (7.08 M allocations: 364.345 MiB, 0.62% compilation time)
elapsed time (ns):  557904420654
gc time (ns):       0
bytes allocated:    382043682
pool allocs:        7070597
non-pool GC allocs: 4410
malloc() calls:     15
realloc() calls:    2
free() calls:       0
minor collections:  0
full collections:   0
160000×80000 SparseMatrixCSC{ComplexF64, Int64} with 320000 stored entries:
⢳⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⢳⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⢳⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⢳⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⢳⡀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⢳⡀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠈⢷⡀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⢷⡀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠈⢷⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠈⢧⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢧⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢧⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢧⠀
⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢧

Second execution is not much different:

julia> @timev A \ B
557.779810 seconds (33 allocations: 19.246 MiB)
elapsed time (ns):  557779810182
gc time (ns):       0
bytes allocated:    20180968
pool allocs:        12
non-pool GC allocs: 13
malloc() calls:     6
realloc() calls:    2
free() calls:       0
minor collections:  0
full collections:   0
160000×80000 SparseMatrixCSC{ComplexF64, Int64} with 320000 stored entries:
⢳⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⢳⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⢳⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⢳⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⢳⡀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⢳⡀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠈⢷⡀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⢷⡀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠈⢷⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠈⢧⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢧⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢧⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢧⠀
⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢧

On the other hand, here is the timing of A \ B for the two matrices constructed by the same code, but on Julia 1.7. First execution takes 3.2 seconds:

julia> VERSION
v"1.7.2-pre.0"

julia> @timev A \ B
  3.199857 seconds (6.39 M allocations: 349.351 MiB, 2.93% gc time, 99.37% compilation time)
elapsed time (ns): 3199857054
gc time (ns):      93635079
bytes allocated:   366321440
pool allocs:       6385404
non-pool GC allocs:4326
malloc() calls:    7
realloc() calls:   2
GC pauses:         2
160000×80000 SparseMatrixCSC{ComplexF64, Int64} with 320000 stored entries:
⢳⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⢳⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⢳⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⢳⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⢳⡀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⢳⡀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠈⢷⡀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⢷⡀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠈⢷⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠈⢧⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢧⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢧⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢧⠀
⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢧

Second execution is even faster with 0.023 secs:

julia> @timev A \ B
  0.023254 seconds (33 allocations: 19.246 MiB)
elapsed time (ns): 23254352
bytes allocated:   20180968
pool allocs:       12
non-pool GC allocs:13
malloc() calls:    6
realloc() calls:   2
160000×80000 SparseMatrixCSC{ComplexF64, Int64} with 320000 stored entries:
⢳⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⢳⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⢳⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⢳⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⢳⡀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⢳⡀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠈⢷⡀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⢷⡀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠈⢷⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠈⢧⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢧⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢧⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢧⠀
⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢧

Is this a known issue? Otherwise, please advise what I should do to further identify the cause of the regression. The two matrices A and B are constructed by my big packages, so it is a bit difficult to generate a MWE demonstrating the regression without further guidance.

Here is the versioninfo() output for Julia 1.8.

julia> versioninfo()
Julia Version 1.8.3-pre.2
Commit 97ccb970c0 (2022-10-28 06:19 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin21.6.0)
  CPU: 8 × Intel(R) Core(TM) i7-6700HQ CPU @ 2.60GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-13.0.1 (ORCJIT, skylake)
  Threads: 1 on 4 virtual cores
Environment:
  JULIA_PKG_DEVDIR = /Users/wsshin/Programming
  JULIA_EDITOR = code
  JULIA_NUM_THREADS = 

[rayegun]

Can you serialize the matrices for me? And upload those files somewhere?

Or some method to generate something with a similar structure.

[SobhanMP]

using SparseArrays, LinearAlgebra

a = spzeros(ComplexF64, 100_000, 100_000) + I
b = spzeros(ComplexF64, 100_000, 50_000)
a \ b

no need to go far

[KristofferC]

In 1.7 the dispatch goes to:

https://github.com/JuliaLang/julia/blob/b06a24a7eed4464e50f4ee594c3a303ad543c42b/stdlib/LinearAlgebra/src/diagonal.jl#L595-L596

which dispatches to

SparseArrays.jl/src/linalg.jl

Lines 892 to 907 in 710fd7b

	function ldiv!(D::Diagonal{T}, A::AbstractSparseMatrixCSC{T}) where {T}
	# require_one_based_indexing(A)
	if size(A, 1) != length(D.diag)
	throw(DimensionMismatch("diagonal matrix is $(length(D.diag)) by $(length(D.diag)) but right hand side has $(size(A, 1)) rows"))
	end
	nonz = nonzeros(A)
	Arowval = rowvals(A)
	b = D.diag
	for i=1:length(b)
	iszero(b[i]) && throw(SingularException(i))
	end
	@inbounds for col in 1:size(A, 2), p in nzrange(A, col)
	nonz[p] = b[Arowval[p]] \ nonz[p]
	end
	A
	end

In 1.8 it first goes to

https://github.com/JuliaLang/julia/blob/0434deb161e17103eabdd7d82b17a1cd6b410572/stdlib/LinearAlgebra/src/diagonal.jl#L422-L423

but the first argument there does not evaluate to a diagonal:

2|julia> _init(\, D, B).(B) |> typeof
SparseMatrixCSC{ComplexF64, Int64}

so we dispatch to

https://github.com/JuliaLang/julia/blob/0434deb161e17103eabdd7d82b17a1cd6b410572/stdlib/LinearAlgebra/src/diagonal.jl#L426-L440

which is not aware of sparsity.

[dkarrasch]

How "embarrassing"... We do have *(S::Sparse, D::Diagonal), *(D, S), and /(S, D) overloads (and corresponding l[mul/div]! and r[mul/div]!, but no \(D, S).

[wsshin]

Is there a way to test the above merge myself? I would like to confirm that the merge fixes my issue.

I have tested merges for regular packages, but I'm not sure how to do it for standard libraries like SparseArrays.

[wsshin]

For example, if I build Julia from HEAD of the master branch, does it automatically retrieve HEAD of the master branch of all the standard libraries?

[rayegun]

No. What you can do is use Julia 1.9-DEV (or master) and then clone this repository. Change the UUID to something else (just change one character) in Project.toml and then if you activate that Project.toml using SparseArrays will use that "local" version, rather than the built-in version.

See this post: https://discourse.julialang.org/t/testing-changes-to-stdlib-without-remaking-the-whole-of-julia/11998/9

[wsshin]

Can this be backported to Julia 1.8? I thought v1.8.3 would be the final Julia 1.8, but I just saw that there is going to be v1.8.4: https://discourse.julialang.org/t/julia-1-8-4-testing-period/91755. I am stuck to Julia 1.7 because of this issue.

[KristofferC]

I think it is backported

[wsshin]

I think we have this issue again in Julia 1.9. Executing @SobhanMP's example in Julia 1.9:

julia> VERSION
v"1.9.0-beta3.1"

julia> @time a \ b;
185.706327 seconds (29 allocations: 7.390 MiB)

On the other hand, in Julia 1.8:

julia> VERSION
v"1.8.5"

julia> @time a \ b;
  0.005846 seconds (29 allocations: 7.390 MiB)

For both cases @time is performed for the second time a \ b is executed.

Can we create a unit test that checks this regression, so that the same bug does not keep reintroduced?

[dkarrasch]

I think we have this issue again in Julia 1.9.

No, the fix hasn't been backported to v1.9 yet, but it is in the pipeline: #300