Update Metis.graph to support only one triangle of a sparse Hermitian matrix

[amontoison]

We could add a keyword argument view in the following routine:

Metis.graph(G::SparseMatrixCSC; check_hermitian=true)

such that the new signature is

Metis.graph(G::SparseMatrixCSC; check_hermitian=true; view='F')

With this modification, we can handle the cases where only the lower or upper triangle part of G is stored.

function graph(G::SparseMatrixCSC; weights::Bool=false, check_hermitian::Bool=true; view::Char='F')
     (view == 'L') || (view == 'U') || (view == 'F') || throw(ArgumentError("The view $view is not supported."))
    if check_hermitian && view == 'F'
        ishermitian(G) || throw(ArgumentError("matrix must be Hermitian"))
    end
    if view == 'L'
      ...
    elseif view == 'U'
      ...
    else
      N = size(G, 1)
      xadj = Vector{idx_t}(undef, N+1)
      xadj[1] = 1
      adjncy = Vector{idx_t}(undef, nnz(G))
      vwgt = C_NULL # TODO: Vertex weights could be passed as input argument
      adjwgt = weights ? Vector{idx_t}(undef, nnz(G)) : C_NULL
      adjncy_i = 0
      @inbounds for j in 1:N
          n_rows = 0
          for k in G.colptr[j] : (G.colptr[j+1] - 1)
              i = G.rowval[k]
              i == j && continue # skip self edges
              n_rows += 1
              adjncy_i += 1
              adjncy[adjncy_i] = i
              if weights
                  adjwgt[adjncy_i] = G.nzval[k]
              end
          end
          xadj[j+1] = xadj[j] + n_rows
      end
      resize!(adjncy, adjncy_i)
      weights && resize!(adjwgt, adjncy_i)
      return Graph(idx_t(N), xadj, adjncy, vwgt, adjwgt)
  end
end

cc @frapac @sshin23

[fredrikekre]

Metis still needs the full graph right? Why not use the existing Hermitian wrapper from LinearAlgebra instead of the view argument?

[amontoison]

I checked a little bit and it seems that Metis always need the full graph.
The issue with the wrapper Hermitian is that both triangles can be also stored so we must be careful if we rely on A.uplo.
We can use the wrapper, but we must verify and discard the values that potentially belong to the other triangle.

[fredrikekre]

What do you mean? If someone is passing Hermitian we should only consider the one half (even if the other one is stored, it should be ignored).

[amontoison]

Yes, you answered before that I modified my last message.
After we can also add a dispatch like that:

Metis.graph(G::Hermitian{<:All, <:SparseMatrixCSC}) = Metis.graph(G.data; check_hermitian=false; view=G.uplo)

[fredrikekre]

What I mean is, what would view bring that you couldn't do with Hermitian?

[amontoison]

We could do everything with Hermitian.
In our optimization solver, we just don't wrap the matrix with Hermitian, so it was great for us to also have a Julia function without the wrapper.

[mipals]

I am having some similar issues. How did you end up solving it @amontoison ?

[amontoison]

I still provide both triangles to METIS. 😕

[mipals]

Okay, I guess I have to do the same then. Not the most efficient for large problems 😬

[mipals]

Okay having looked at the code and the Metis documentation I think the issue is not as bad as I initially thought. Yes, we need to store each edge twice (to represent the full sparsity pattern), but the weights are not needed. I guess this is what Alexis hinted at earlier with the specialized graph assembly for when only tril or triu is present.

[amontoison]

It's what we do in the HSL linear solvers like MA57, the user only provides one triangle and after AMD or METIS is called to compute a fill-in reduction permutation.
I should check the Fortran code to see how they build the graph for METIS.

[mipals]

Yeah, that is pretty standard for sparse direct solvers. You also only provide one of the triangles in Pardiso and faer-rs.

[mipals]

My rough solution so far is to get the sparsity pattern of the missing triangle by converting the original triangle using a "CSR to CSC" procedure without copying the non-zeros (the transpose of a CSC is simply a CSR with colptr and rowval as rowptr and colval if that makes sense). The final step is then to merge the two triangular sparsity patterns, which can be done fairly straightforward using vcat because of the triangular structure .

[mipals]

A rough implementation of the above can be seen in #54 . Any diagonal entries will be duplicated and subsequently ignored when creating the Graph. This could probably be done more elegantly.

[fredrikekre]

Closed by #54.