Should we check if row indices is sorted in every column when importing sparse matrix?

[songjhaha]

In docs of SparseArrays.jl

If you have data in CSC format from a different application or library, and wish to import it in Julia, make sure that you use 1-based indexing. The row indices in every column need to be sorted, and if they are not, the matrix will display incorrectly. If your SparseMatrixCSC object contains unsorted row indices, one quick way to sort them is by doing a double transpose. Since the transpose operation is lazy, make a copy to materialize each transpose.

And i get some .mat data from https://github.com/qpsolvers/maros_meszaros_qpbenchmark/blob/main/data/QFORPLAN.mat which contains sparse matrix with unsorted row indices. Then the getindex method would crash because Julia remove inbounds check.

We could check the Ir and Jc to make sure row indices is sorted, If it's unsorted, use double transpose to collect the correct result?

It's a feature request and I would make an PR if we agree it should be done in MAT.jl instead of SparseArrays.jl

[matthijscox]

I suppose I'm okay with doing that here in MAT.jl if SparseArrays.jl doesn't provide any checks/functionality for this. I admit I have not touched the sparse array code yet since becoming maintainer.

Just to be clear; this file currently loads without error in MAT.jl, but the resulting matrix is corrupt/incorrect when you start indexing?

I just tried and reading works indeed. I see two sparse matrices (A and P), and I seem capable of indexing them?

julia> vars = matread("QFORPLAN.mat")
Dict{String, Any} with 8 entries:
  "A" => sparse([2, 162, 2, 163, 3, 164, 3, 165, 4, 166  …  47, 578, 46, 579, 48, 580, 45, 581, 44, 582], [1, 1, 2, 2, …
  "m" => 582.0
  "P" => sparse([1, 2, 62, 63, 64, 65, 66, 67, 68, 69  …  74, 75, 76, 77, 78, 79, 80, 81, 82, 83], [1, 1, 1, 1, 1, 1, 1…
  "r" => 0.0
  "l" => [7.392e6; -3.55271e-15; … ; 0.0; 0.0;;]
  "q" => [0.02466; 0.0; … ; -0.0012632; -0.0010105;;]
  "u" => [7.392e6; -3.55271e-15; … ; 1.0e20; 1.0e20;;]
  "n" => 421.0

julia> vars["P"][400,400]
0.0

[songjhaha]

It would corrupt when slice vars["A"][1:end,:]

[matthijscox]

Still not sure what you mean, indexing like that works fine for me:

julia> vars["A"][1:end,:] == vars["A"]
true

Converting to non-sparse matrices/vectors also works fine:

julia> Matrix(vars["A"][1:end,:]) == Matrix(vars["A"])
true

julia> Vector(vars["A"][3,:]) == Matrix(vars["A"])[3,:]
true

Bounds checks also work:

julia> vars["A"][583,:]
ERROR: BoundsError: attempt to access 582×421 SparseArrays.SparseMatrixCSC{Float64, Int64} with 4984 stored entries at index [583, 1:421]

[songjhaha]

sorry, here is the MWE:

> julia --check-bounds=yes

using MAT
vars = matread("QFORPLAN.mat")
vars["A"][vec(vars["l"] .== vars["u"]),:]

result:

ERROR: BoundsError: attempt to access 3247-element Vector{Int64} at index [3248]
Stacktrace:
 [1] setindex!
   @ .\array.jl:1021 [inlined]
 [2] getindex_I_sorted_linear(A::SparseArrays.SparseMatrixCSC{Float64, Int64}, I::Vector{Int64}, J::UnitRange{Int64})
   @ SparseArrays C:\Users\Public\TongYuan\julia-1.10.10\share\julia\stdlib\v1.10\SparseArrays\src\sparsematrix.jl:2905
 [3] getindex_I_sorted(A::SparseArrays.SparseMatrixCSC{Float64, Int64}, I::Vector{Int64}, J::UnitRange{Int64})
   @ SparseArrays C:\Users\Public\TongYuan\julia-1.10.10\share\julia\stdlib\v1.10\SparseArrays\src\sparsematrix.jl:2792
 [4] getindex(A::SparseArrays.SparseMatrixCSC{Float64, Int64}, I::Vector{Int64}, J::UnitRange{Int64})
   @ SparseArrays C:\Users\Public\TongYuan\julia-1.10.10\share\julia\stdlib\v1.10\SparseArrays\src\sparsematrix.jl:3060
 [5] getindex
   @ C:\Users\Public\TongYuan\julia-1.10.10\share\julia\stdlib\v1.10\SparseArrays\src\sparsematrix.jl:3118 [inlined]
 [6] getindex(A::SparseArrays.SparseMatrixCSC{Float64, Int64}, i::BitVector, ::Colon)
   @ SparseArrays C:\Users\Public\TongYuan\julia-1.10.10\share\julia\stdlib\v1.10\SparseArrays\src\sparsematrix.jl:2688
 [7] top-level scope
   @ REPL[4]:1

[matthijscox]

Alright, I was able to reproduce this.

The double copy transpose did allow indexing:

idx = findall(vec(vars["l"] .== vars["u"]))
A = vars["A"]
At = copy(copy(A')')
At[idx,:] # this works

Not sure if we can avoid the double copying somehow. Also note that we read SparseMatrixCSC in both MAT_v5.jl and MAT_HDF5.jl.

Feel free to open a PR to try to fix this.

[songjhaha]

I could try some sorted method like https://github.com/scipy/scipy/blob/8a4633fa0e01d62e9ccdd06ebe5bb30551cfa056/scipy/sparse/sparsetools/csr.h#L368

And compare the performance with double transpose method