Adding support for sparse Hermitian matrices.

[mipals]

Initial idea for computing the full sparsity pattern of a sparse Hermitian matrix by using only its lower or upper triangular part.

[fredrikekre]

This algorithm actually exist in SparseArrays already for converting Hermitian(SparseMatrix) to SparseMatrix, see https://github.com/JuliaSparse/SparseArrays.jl/blob/313a04f4a78bbc534f89b6b4d9c598453e2af17c/src/sparseconvert.jl#L124-L173

Here is a modified version of that which ignores the diagonal

const HermOrSymCSC{Tv, Ti} = Union{
    Hermitian{Tv, SparseMatrixCSC{Tv, Ti}}, Symmetric{Tv, SparseMatrixCSC{Tv, Ti}},
}

"""
    Metis.graph(G::Union{Hermitian, Symmetric}; weights=false)

Construct the 1-based CSR representation of the `Hermitian` or `Symmetric` wrapped sparse
matrix `G`.
Weights are not currently supported for this method so passing `weights=true` will throw an
error.
"""
function graph(H::HermOrSymCSC; weights::Bool=false)
    weights && throw(ArgumentError("weights not supported yet"))
    # Extract data
    A = H.data
    upper = H.uplo == 'U'
    rowval = rowvals(A)
    m, n = size(A)
    @assert m == n
    # New colptr for the full matrix
    newcolptr = Vector{idx_t}(undef, n + 1)
    newcolptr[1] = 1
    # SparseArrays.nzrange for the upper/lower part excluding the diagonal
    nzrng = if upper
        (A, col) -> SparseArrays.nzrangeup(A, col, #=exclude diagonal=# true)
    else
        (A, col) -> SparseArrays.nzrangelo(A, col, #=exclude diagonal=# true)
    end
    # If the upper part is stored we loop forward, otherwise backwards
    colrange = upper ? (1:1:n) : (n:-1:1)
    @inbounds for col in colrange
        r = nzrng(A, col)
        # Number of entries in the stored part of this column, excluding the diagonal entry
        newcolptr[col + 1] = length(r)
        # Increment columnptr corresponding to the stored rows
        for k in r
            row = rowval[k]
            @assert upper ? row < col : row > col
            @assert row != col # Diagonal entries should not be here
            newcolptr[row + 1] += 1
        end
    end
    # Accumulate the colptr and allocate new rowval
    cumsum!(newcolptr, newcolptr)
    nz = newcolptr[n + 1] - 1
    newrowval = Vector{idx_t}(undef, nz)
    # Populate the rowvals
    @inbounds for col = 1:n
        newk = newcolptr[col]
        for k in nzrng(A, col)
            row = rowval[k]
            @assert col != row
            newrowval[newk] = row
            newk += 1
            ni = newcolptr[row]
            newrowval[ni] = col
            newcolptr[row] = ni + 1
        end
        newcolptr[col] = newk
    end
    # Shuffle back the colptrs
    @inbounds for j = n:-1:1
        newcolptr[j+1] = newcolptr[j]
    end
    newcolptr[1] = 1
    # Return Graph
    N = n
    xadj = newcolptr
    adjncy = newrowval
    vwgt = C_NULL
    adjwgt = C_NULL
    return Graph(idx_t(N), xadj, adjncy, vwgt, adjwgt)
end

How does this compare to your version? This one only allocates the vectors once at least.

[mipals]

I did actually look in SparseArrays for something, but with no luck. Thats one of the reason I wrote on Slack - In the hopes that someone knew something 😄

This seems seem work, and it is definitely more efficient that what I had as it does not allocate things twice and ignores the diagonal directly 👍

[fredrikekre]

Thanks, available in version 1.5.0 (JuliaRegistries/General#115506).

[amontoison]

Thanks @mipals for solving #53!

[mipals]

In the end I think @fredrikekre ended up doing the work. I think I just pushed him in the right direction. So the thanks should really go to him 🥳