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. With the max_diagonal_product strategy it also emits row and column scaling coefficients \(D_r, D_c\) 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\) (permutation only)

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.

For max_diagonal_product, MC64 chooses \(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\). The max_diagonal_sum strategy computes the permutation only; scaling coefficients are computed for max_diagonal_product.

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<ValueType, IndexType>::generate returns a gko::Composition of two gko::matrix::ScaledPermutation objects. The first holds the row scaling and the row permutation; the second holds the column scaling (with the identity as its permutation). Csr::scale_permute(row, col) applies both to a matrix in one call.

using mc64 = gko::experimental::reorder::Mc64<double, int>;
using scaled_perm = gko::matrix::ScaledPermutation<double, int>;
auto mc64_factory =
    mc64::build()
        .with_strategy(
            gko::experimental::reorder::mc64_strategy::max_diagonal_product)
        .on(exec);
auto result = mc64_factory->generate(system_matrix);

auto row_perm = gko::as<scaled_perm>(result->get_operators()[0]);
auto col_perm = gko::as<scaled_perm>(result->get_operators()[1]);
auto matched = system_matrix->scale_permute(row_perm, col_perm);

Parameter	Type	Purpose
`strategy`	`mc64_strategy`	Goal of the matching: `max_diagonal_product` (default) or `max_diagonal_sum`.
`tolerance`	`remove_complex<ValueType>`	Tolerance below which a weight is considered zero. Defaults to `50 * epsilon` of the real value type.

The matched matrix can then be reordered with AMD or nested dissection and factorised; see Reordering and permutations for the full recipe.

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.