Bug / Design Mismatch: Sparse LU Fill-Reduction Permutations

[micycle1]

While preparing a PR to add SYMRCM as a fill-reduction option for sparse CSC LU, I found what looks like a design/implementation mismatch in the current LU fill-permutation path, and I wanted to check the intended approach before proceeding.

What I tested

I added a symmetric Reverse Cuthill-McKee ordering (SYMRCM) and tried using it via:

LinearSolverFactory_DSCC.lu(FillReducing.SYMRCM)

On a square matrix with a structurally symmetric sparsity pattern:

• IDENTITY works
• SYMRCM gives incorrect results

What seems to be going on

The LU path appears to mix two different models for fill-reduction permutations.

1. ApplyFillReductionPermutation_DSCC

This class is documented as applying P*A*Q.

But for the nonsymmetric case it only physically applies a row permutation:

CommonOps_DSCC.permuteRowInv(pinv, A, Aperm);

and keeps the column permutation separate.

It also requires a row permutation to exist; if I try a column-only fill permutation, LU fails immediately with:

RuntimeException: No row permutation matrix

2. LuUpLooking_DSCC.performLU()

Inside LU, the factorization loop uses:

int[] q = applyReduce.getArrayP();
int col = q != null ? q[k] : k;

So it is taking getArrayP() (row permutation) and using it as the column ordering for LU. For SYMRCM, where prow == pcol, that happens to mask the issue, but semantically it looks wrong; I would have expected getArrayQ() here.

3. LinearSolverLu_DSCC.solve()

The solve phase does:

• inverse numeric row pivot (pinv)
• forward/back substitution
• apply one reduce permutation back to the solution

i.e. it clearly accounts for:

• numeric LU row pivoting
• one fill-reduction permutation on the solution side

But it does not appear to apply any fill-reduction row permutation to the RHS before solving.

So if LU is factorizing something equivalent to:

$$ P_f A Q_f $$

then solve should also account for P_f b, but that does not seem to happen.

Why I think this matters

This suggests the current LU fill-reduction path is internally inconsistent for non-identity permutations:

• it requires a row permutation,
• but solve only visibly handles the column side of fill reduction,
• and LU appears to use P where Q would be expected.

That would explain why:

• IDENTITY works,
• but nontrivial fill reduction like SYMRCM does not.

Question

Before I proceed with a PR, I wanted to ask what the intended model for sparse LU fill reduction is supposed to be:

Option A

LU should support only a column fill-reduction permutation, with row permutations handled only by numeric pivoting.

In that case, the LU path should probably:

• stop requiring a row permutation in ApplyFillReductionPermutation_DSCC
• use getArrayQ() in LU
• keep the existing solution-side unpermutation

Option B

LU should support a full static fill permutation P*A*Q.

In that case, solve would also need to apply the fill row permutation to the RHS, and LU should consistently distinguish P from Q.

---

The core question is: for sparse LU in EJML, is fill reduction intended to be column-only, or full P*A*Q?

[lessthanoptimal]

Thanks for looking into this. I'm going to need to knock off some rust and dig into this code to refresh my memory. Have you looked at the unit tests? In theory, they should go over these cases. This is also very well tread territory, so any PR should cite a source for the algorithm. Just added this to my todo list and will hopefully get back to you in a few days.