KroneckerStencilMatrix: grouped factors, ComposedKroneckerStencilMatrix, KroneckerSumSolver

[spossann]

Summary

Improves KroneckerStencilMatrix and KroneckerLinearSolver in feectools/linalg/kron.py.

Factors on groups of axes. A factor of a KroneckerStencilMatrix may now act on several consecutive axes, e.g. a 2d x 1d product on a 3d space:

M = KroneckerStencilMatrix(V, W, A_xy, A_z)   # M.axes == ((0, 1), (2,))

• New axes property. ndim is now the number of axes of the domain; before, it was the number of factors, which is the same for 1d factors.
• Stricter constructor checks: the factors must be StencilMatrix, the codomain npts are checked as well, and the factors must cover exactly V.ndim axes.
• dot, tostencil and __getitem__ work with grouped factors. dot and tostencil now use each factor's own pads and row offsets, so a factor's pads may be smaller than the pads of the space.

Matrix product A @ B → ComposedKroneckerStencilMatrix. For two Kronecker matrices with the same axis groups, A @ B returns a new ComposedKroneckerStencilMatrix. This mirrors LinearOperator @ LinearOperator → ComposedLinearOperator, and the new class subclasses ComposedLinearOperator:

• multiplicands are the operands. Chains are flattened on both sides, so (A @ B) @ C and A @ (B @ C) both have (A, B, C). dot applies them from right to left, as in ComposedLinearOperator.
• mats holds the exact Kronecker factors C_k = A_k @ B_k of the product. They are computed once at construction in scipy sparse format and stored on process-local spaces with wider pads (p_A + p_B). If the operands are distributed, the rows of B_k owned by other processes are gathered, so construction is collective.
• tosparse and toarray use the exact C_k, which is also correct in parallel. A KroneckerLinearSolver for the product can be built from mats.
• Scaling, copy, transpose and further products with Kronecker matrices keep the type. Any other operand (another operator type, other axis groups) gives a plain ComposedLinearOperator.
• Naming: mats always means the Kronecker factors (A₁ ⊗ A₂ ⊗ …), and multiplicands the operands of the matrix product (A · B · …).

KroneckerStencilMatrix itself takes no extra argument for products. Its factor pads must always fit into the domain pads, so dot and tostencil always work.

multiplicants → multiplicands. ComposedLinearOperator.multiplicants is renamed to multiplicands, the correct spelling. multiplicants remains as a deprecated alias that emits a DeprecationWarning. All uses in feectools (api/fem_common.py, api/fem_bilinear_form.py) are updated.

KroneckerLinearSolver / kronecker_solve. New optional argument factor_ndims (default: all 1, so existing behaviour is unchanged). A solver of a factor with several axes receives the vectors flattened over these axes in C order, as in StencilMatrix.tosparse. The axes of such a factor must not be distributed across processes; otherwise NotImplementedError is raised. 1d factors are solved in parallel as before.

KroneckerSumSolver (fast diagonalization). A new solver for sums of Kronecker products,

A = Σ_d M_1 ⊗ … ⊗ S_d ⊗ … ⊗ M_n + σ M_1 ⊗ … ⊗ M_n

with symmetric 1d matrices S_d and symmetric positive definite 1d matrices M_d, e.g. a Laplacian on a tensor-product grid. Such a sum is not a single Kronecker product, so KroneckerLinearSolver can't invert it.

• In each direction, the generalized eigenproblem S_d U_d = M_d U_d Λ_d is solved (with U_dᵀ M_d U_d = I). Then A⁻¹ = U Λ⁻¹ Uᵀ, with U = U_1 ⊗ … ⊗ U_n and the diagonal Λ = Λ_1 ⊕ … ⊕ Λ_n + σ.
• Uᵀ and U are applied with KroneckerLinearSolver, through a small LinearSolver that multiplies by a dense matrix, as in fft.py. So it also works along distributed axes. Λ⁻¹ is applied locally.
• Vanishing eigenvalues are skipped (pseudo-inverse), e.g. the constants of a periodic Laplacian with σ = 0.
• A direction without a stiffness term is given as None.

It is used by the new StiffnessPreconditioner in struphy (see the paired PR).

DirectionalDerivativeOperator gets the public properties diffdir, negative and transposed.

Other

• Fixed KroneckerStencilMatrix.__imul__, which assigned into a tuple and raised TypeError.
• Docstrings (numpydoc) and type annotations for both classes and kronecker_solve.
• Version bump to 0.7.0 (devel-tiny, which is merged into this branch, is already at 0.6.0).
• compile_psydac.mk: pyccel compile without -v, for a shorter compile output.

Tests

New tests in linalg/tests/test_kron_stencil_matrix.py compare against dense matrices for the groupings (1,1,1), (2,1), (1,2) and (3,). They cover:

• dot, tosparse, __getitem__, scaling and tostencil for grouped factors;
• @: type, values, multiplicands, chains on both sides, scaling, copy, transpose, the fallback to ComposedLinearOperator, and the ValueError case;
• the solver with factor_ndims, including a solver for a product;
• the deprecated multiplicants alias;
• KroneckerSumSolver against dense solves, regular and singular (pseudo-inverse), with and without a direction that has no stiffness term.

Run locally:

• serial: 19 passed;
• mpirun -n 2 and -n 4: all Kronecker and test_fft MPI tests pass, including distributed axes for KroneckerSumSolver;
• full serial feectools suite, with feectools DeprecationWarnings turned into errors: 9101 passed.

The four collection errors in ddm/tests/test_cart_2d.py, test_cart_3d.py, linalg/tests/test_block.py and test_toarray.py (unregistered parallel marker) also occur without this change.

Transposes are tested in serial only: in parallel, the process-local factors don't hold the rows of other processes. This was already the case before this PR.

Paired struphy PR (runs the struphy tests against this branch): struphy-hub/struphy#731

🤖 Generated with Claude Code

[Copilot]

🟡 Changes recommended

Parallel products can discard valid periodic entries, and incompatible process-local factor layouts are not rejected.

2 open findings

• Validate factor codomain row layout matches W · New
• Canonicalize COO duplicates before MPI entry removal · New

What changed in this PR

Extends Kronecker linear algebra with grouped factors, composed products, and fast diagonalization.

Changes:

• Adds grouped-axis Kronecker matrices and composed products.
• Adds grouped-factor and Kronecker-sum solvers.
• Renames multiplicants and exposes derivative metadata.

File	Description
pyproject.toml	Bumps version to 0.6.0.
feectools/​linalg/​kron.py	Implements the new Kronecker functionality.
feectools/​linalg/​basic.py	Renames composed-operator multiplicands.
feectools/​feec/​derivatives.py	Adds derivative properties.
feectools/​api/​fem_common.py	Uses the renamed property.
feectools/​api/​fem_bilinear_form.py	Uses the renamed property.
feectools/​linalg/​tests/​test_kron_stencil_matrix.py	Adds extensive serial and MPI coverage.

🧠 Review effort: Balanced

---

Give feedback about Copilot approvals in this survey to enter a drawing for a $150 gift card.

[spossann]

@max-models this is ready for review!