Sparse matrices: Catastrophic fill-in when solving transpose

[eschnett]

I am solving a linear system where the LHS is a diagonal matrix and the RHS is the transpose of a sparse matrix. The resulting sparse matrix has all its zeros filled in. I show an example where I first create a diagonal LHS and a sparse RHS, and the solve with both the sparse matrix and its transpose. Note the with 9 stored entries in the output; this should not be.

Brief testing indicates that this is also related to rational numbers of big integers.

julia> diag = Diagonal([big(1)//1,1,1])
3×3 Diagonal{Rational{BigInt}, Vector{Rational{BigInt}}}:
 1  ⋅  ⋅
 ⋅  1  ⋅
 ⋅  ⋅  1

julia> sp = sparse(Diagonal([1,1,1]))
3×3 SparseMatrixCSC{Int64, Int64} with 3 stored entries:
 1  ⋅  ⋅
 ⋅  1  ⋅
 ⋅  ⋅  1

julia> diag \ sp
3×3 SparseMatrixCSC{Rational{BigInt}, Int64} with 3 stored entries:
 1  ⋅  ⋅
 ⋅  1  ⋅
 ⋅  ⋅  1

julia> diag \ sp'
3×3 SparseMatrixCSC{Rational{BigInt}, Int64} with 9 stored entries:
 1  0  0
 0  1  0
 0  0  1

julia> versioninfo()
Julia Version 1.10.0-beta2
Commit a468aa198d0 (2023-08-17 06:27 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: macOS (x86_64-apple-darwin22.4.0)
  CPU: 16 × Intel(R) Core(TM) i9-9980HK CPU @ 2.40GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, skylake)
  Threads: 1 on 16 virtual cores

[jishnub]

Seems related to JuliaSparse/SparseArrays.jl#389

julia> B = spzeros(3,3); B[1,1] = 0.0; B
3×3 SparseMatrixCSC{Float64, Int64} with 0 stored entries:
  ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅ 

julia> B = spzeros(BigFloat,3,3); B[1,1] = big(0.0); B
3×3 SparseMatrixCSC{BigFloat, Int64} with 1 stored entry:
 0.0   ⋅    ⋅ 
  ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅ 

Although these methods perhaps should be specialized in SparseArrays, as these may be optimized considerably by using the structure.