ralna / spral

Sparse Parallel Robust Algorithms Library
https://ralna.github.io/spral/
Other
104 stars 27 forks source link

Avoid hardcoding METIS 5 options #134

Closed jfowkes closed 1 year ago

jfowkes commented 1 year ago

Resolves #133 in a fully backward compatible way by dropping the METIS_OPTION_NUMBERING=1 option (Fortran-style 1-based indexing) and instead using the METIS 5 default (C-style 0-based indexing) and adding 1 to the perm and invp arrays returned by METIS 5 (the fill-reducing permutation and inverse permutation respectively) as well as subtracting 1 from the ptr and row input arrays (matrix sparsity pattern in CSR format).

Note that there is virtually no performance loss with this change as the perm and invp arrays returned by METIS 5 are already copied in the wrapper (this is so that their return types are correct). The same is true of the input arrays ptr and row that encode the matrix sparsity pattern in CSR format which are copied so that the matrix is stored fully (not merely its lower-triangular part) as is expected by METIS 5.

jfowkes commented 1 year ago

@mjacobse I'd be interested in your thoughts on this proposed fix, does it make sense to you to do it this way?

mjacobse commented 1 year ago

Makes sense to me 👍

One thing I would like to suggest is to not do the initial conversion from Fortran to C indexing in the implementations of the half_to_full_drop_diag interface, but instead after the calls to it in metis_order32/64. That would reduce the places in the code where this is done from 4 to 2. It would also make half_to_full_drop_diag do only the one thing that the name suggests, and not an additional somewhat unexpected index shift. And it would bring the forward and backward index-shift into the same function, right next to the Metis call for which it is actually done. Both places should even fit on screen at the same time, so it would be much easier to directly see and reason about while reading the code.

I'm not quite sure what the change from dimension(*) to dimension(:) is for in half_to_full_drop_diag? Could this also be dropped if moving the index shift outside of it?

Alternatively, if the index shift is to remain within half_to_full_drop_diag, then I think the documentation and perhaps also names should be updated accordingly. It could then also be considered to do the index shift directly while building the row2 and ptr2 arrays instead of afterwards.

jfowkes commented 1 year ago

Thank you, yes excellent point agreed, I will move the initial conversion over to metis_order32/64 instead of doing it in half_to_full_drop_diag which is messy.

The change from dimension(*) to dimension(:) in half_to_full_drop_diag was only there so that I could subtract one from the arrays without specifying extents so yes this can be dropped for consistency with the METIS 4 wrapper.