MC64 (matching and scaling)#

MC64 is a matching-and-scaling preprocessing step for sparse direct factorisation. Unlike AMD or nested dissection (which permute rows and columns to reduce fill-in but leave numerical values untouched), MC64 picks a row permutation that maximises a metric over the resulting diagonal absolute values and simultaneously computes a pair of diagonal scalings \(D_r, D_c\) (for max_diagonal_product) that bring every diagonal entry of the transformed matrix to unit magnitude.

Phrased as an optimisation, MC64 solves a bipartite-matching problem on the matrix graph. The two scoring strategies Ginkgo supports correspond to different per-edge weight functions \(w_{ij}\), and in both cases the matching minimises the reduced cost \(m_i - w_{i,\sigma(i)}\), where \(m_i = \max_j w_{ij}\) is the row-maximum:

Strategy	Per-edge weight \(w_{ij}\)	Matching objective	Equivalent diagonal goal
`max_diagonal_product` (default)	\(w_{ij} = \log_2 \lvert a_{ij} \rvert\)	\(\min_{\sigma} \sum_{i} \bigl( m_i - w_{i,\sigma(i)} \bigr)\)	maximise \(\prod_{i} \lvert a_{i,\sigma(i)} \rvert\)
`max_diagonal_sum`	\(w_{ij} = \lvert a_{ij} \rvert\)	\(\min_{\sigma} \sum_{i} \bigl( m_i - w_{i,\sigma(i)} \bigr)\)	maximise \(\sum_{i} \lvert a_{i,\sigma(i)} \rvert\)

Subtracting the row maximum before matching keeps every cost non-negative, as required by the shortest-augmenting-path matching the implementation uses. The log-space transform on the product strategy also avoids the underflow that would otherwise plague a literal \(\prod \lvert a_{ij} \rvert\) over thousands of subunit values.

With max_diagonal_product, MC64 also computes \(D_r\) and \(D_c\) so the permuted-and-scaled matrix \(\hat A = P D_r\, A\, D_c\) has \(\lvert \hat a_{ii} \rvert = 1\) and \(\lvert \hat a_{ij} \rvert \le 1\); max_diagonal_sum computes no scaling coefficients.

When to use MC64#

Indefinite or strongly non-symmetric matrices headed for LU factorisation. Ginkgo’s Lu factorises without pivoting; pushing large entries onto the diagonal can increase the numerical stability of that factorisation.
As the match step of the canonical recipe match → reorder → factorise; see Reordering and permutations.

For SPD systems factorised by Cholesky, MC64 is best skipped: the diagonal is already heavy, and MC64’s row permutation would break the symmetry Cholesky relies on. For Krylov iterative methods, MC64’s scaling is occasionally useful as a diagonal preconditioner, but the matching step adds value only when the matrix has zero or near-zero diagonal entries.

Using MC64#

Mc64::generate returns a gko::Composition<ValueType> of two gko::matrix::ScaledPermutation<ValueType, IndexType> operators: the row permutation with its row scaling, and the column scaling (with an identity permutation). Csr::scale_permute applies both to a matrix.

using ScaledPerm = gko::matrix::ScaledPermutation<double, int>;
auto mc64 = gko::experimental::reorder::Mc64<double, int>::build()
                .with_strategy(
                    gko::experimental::reorder::mc64_strategy::max_diagonal_product)
                .on(exec)
                ->generate(system_matrix);
auto row = gko::as<ScaledPerm>(mc64->get_operators()[0]);
auto col = gko::as<ScaledPerm>(mc64->get_operators()[1]);
auto scaled = system_matrix->scale_permute(row, col);

The factory parameters are:

Parameter	Type	Purpose
`strategy`	`mc64_strategy`	Matching objective: `max_diagonal_product` (default) or `max_diagonal_sum`.
`tolerance`	`remove_complex<ValueType>`	Weights below this value count as zero; default `50 * epsilon`.

Combine the result with a fill-reducing reordering by applying Csr::permute to the scaled matrix before the factorisation.

Implementation notes#

The algorithm runs on host. The matching is solved by an augmenting-path / Hungarian-style algorithm; expect runtime to be modest for matrices up to a few million nonzeros but to dominate the pre-factorisation budget on much larger problems.