Skip to content

KroneckerStencilMatrix: grouped factors, ComposedKroneckerStencilMatrix, KroneckerSumSolver - #100

Merged
spossann merged 7 commits into
devel-tinyfrom
kron-stencil-matmul
Oct 8, 2026
Merged

KroneckerStencilMatrix: grouped factors, ComposedKroneckerStencilMatrix, KroneckerSumSolver#100
spossann merged 7 commits into
devel-tinyfrom
kron-stencil-matmul

Conversation

@spossann

@spossann spossann commented Oct 8, 2026 •
edited
Loading

Copy link
Copy Markdown
Member

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

@spossann @claude
KroneckerStencilMatrix: matmul, factors on groups of axes, docs
6ee8f08 
- KroneckerStencilMatrix factors may act on several consecutive axes
  (e.g. 2d x 1d on a 3d space); new `axes` property, stricter checks
  (StencilMatrix factors, codomain npts, number of axes).
- `A @ B` of two Kronecker matrices with the same axis groups returns a
  KroneckerStencilMatrix with factors C_k = A_k @ B_k (computed in sparse
  format on process-local spaces with wider pads, rows of B gathered
  across processes). Domain/codomain stay the original spaces; copies of
  the operands are kept in `factors` and used by `dot`.
- dot/tostencil use the factor's own pads; tostencil raises ValueError
  if the band does not fit into the domain pads.
- KroneckerLinearSolver and kronecker_solve take `factor_ndims` for
  solvers of factors with several axes (serial along grouped axes).
- Fix __imul__ on the tuple of factors; docstrings and type annotations.
- Bump version to 0.6.0.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@spossann spossann mentioned this pull request Oct 8, 2026
Merged
@spossann @claude
ComposedKroneckerStencilMatrix for products; rename multiplicands
fa6fe90 
- A @ B of Kronecker matrices with the same axis groups now returns a
  ComposedKroneckerStencilMatrix (subclass of ComposedLinearOperator):
  `multiplicands` are the operands (dot goes through them), `mats` the
  exact factors C_k of the product (process-local, wide band; used by
  tosparse and for KroneckerLinearSolver). Scaling, copy, transpose and
  further products keep the type.
- KroneckerStencilMatrix loses the `factors` argument: its factor pads
  always fit into the domain pads, so dot and tostencil always work.
- ComposedLinearOperator.multiplicants -> multiplicands (correct
  spelling); `multiplicants` remains as a deprecated alias.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@spossann spossann changed the title KroneckerStencilMatrix: matmul, factors on groups of axes, docs KroneckerStencilMatrix: grouped factors, ComposedKroneckerStencilMatrix, multiplicands Oct 8, 2026
spossann and others added 2 commits October 8, 2026 12:05
@spossann @claude
Add KroneckerSumSolver (fast diagonalization)
0546a02 
Exact inverse of sum_d M_1 x ... x S_d x ... x M_n + sigma M_1 x ... x M_n
(e.g. a Laplacian on a tensor-product grid): per direction the generalized
eigenproblem S_d U_d = M_d U_d Lambda_d is solved, and A^{-1} = U Lambda^{-1} U^T
with U = U_1 x ... x U_n. U^T and U are applied with KroneckerLinearSolver
(also along distributed axes), Lambda^{-1} locally; vanishing eigenvalues are
skipped (pseudo-inverse). A direction without stiffness term is given as None.

Tests against dense solves, regular and singular, serial and with MPI.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@spossann @claude
DirectionalDerivativeOperator: diffdir, negative and transposed prope…
c774ecf 
…rties

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@spossann spossann changed the title KroneckerStencilMatrix: grouped factors, ComposedKroneckerStencilMatrix, multiplicands KroneckerStencilMatrix: grouped factors, ComposedKroneckerStencilMatrix, KroneckerSumSolver Oct 8, 2026
@spossann
spossann requested a balanced review from Copilot October 8, 2026 11:35
Copilot AI reviewed Oct 8, 2026
View reviewed changes

Copilot AI left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🟡 Changes recommended

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

2 open findings
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.

Comment thread feectools/linalg/kron.py
Comment thread feectools/linalg/kron.py
@spossann @claude
Kronecker matrices: factor row layout and periodic duplicates (review)
e000fb4 
- KroneckerStencilMatrix: the rows owned by a factor must contain the rows
  of the codomain on this process (ValueError otherwise). dot and tostencil
  index the factor rows by global row, so both process-local factors and
  factors owning all rows (e.g. from tokronstencil) work in parallel;
  before, other rows were silently used.
- Products: sum duplicate COO entries of the local factor (periodic factors
  with 2p + 1 > n have two diagonals in the same column) before gathering
  the rows of other processes; before, np.unique dropped one of them in
  parallel.
- Tests for full-row factors, the row check and periodic duplicates (serial
  and MPI).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@spossann
spossann requested a review from max-models October 8, 2026 12:58
@spossann

spossann commented Oct 8, 2026

Copy link
Copy Markdown
Member Author

@max-models this is ready for review!

@spossann
spossann merged commit 5385af9 into devel-tiny Oct 8, 2026
9 checks passed
spossann added a commit to struphy-hub/struphy that referenced this pull request Oct 9, 2026
@spossann @claude
Kronecker preconditioners: StiffnessPreconditioner, refactor of the m…
0432020 
…ass preconditioners; bump feectools (#731)

## Summary

Points the `feectools` submodule to the branch of
struphy-hub/feectools#100, so that the struphy test suite runs against
it before that PR is merged.

struphy-hub/feectools#100 changes `KroneckerStencilMatrix` and
`KroneckerLinearSolver`:
- factors on groups of axes (e.g. 2d x 1d);
- `A @ B` returns a new `ComposedKroneckerStencilMatrix` (operands in
`multiplicands`, exact Kronecker factors of the product in `mats`);
- `ComposedLinearOperator.multiplicants` → `multiplicands`, with a
deprecated alias;
- `KroneckerSumSolver`, the exact inverse of sums of Kronecker products
by fast diagonalization;
- `factor_ndims` in `KroneckerLinearSolver`;
- docstrings and type annotations; feectools 0.7.0.

Existing calls are unchanged, so struphy needs no code changes. struphy
builds these objects in `feec/preconditioner.py`
(`MassMatrixPreconditioner`) and `feec/projectors.py`.

### `multiplicants` → `multiplicands`

- All struphy uses (preconditioner, projectors, multigrid coarsening)
now use `multiplicands`.
- `linop_nbytes` in `feec/memory.py` adds up all child attributes, so a
`ComposedKroneckerStencilMatrix` counts both its operands and the
factors of the product.

### Refactor of `MassMatrixPreconditioner`

- The construction of the Kronecker approximation is split into
module-level helpers in `feec/preconditioner.py`. They cover the
boundary conditions, gathering array weights over MPI, the 1d weight,
the 1d mass matrix and its BCs, the 1d solver, the process-local 1d
factor, the block-diagonal assembly and the operator to invert.
- `MassMatrixDiagonalPreconditioner` now calls these helpers instead of
repeating the same ~200 lines.
- Values and the order of MPI calls are unchanged. The copy into the
process-local 1d stencil matrix is now vectorized, and the
block-diagonal operators are no longer hard-coded to 3 components.
- Docstrings, comments and type annotations.

### `MassMatrixPreconditioner`: new option and fixes

- New option `weight_reduction="midpoint" | "average"` sets how the
weight is reduced to 1d in the directions other than `dim_reduce`. The
default `"midpoint"` gives the same results as before, verified for M0,
M1, M2, Mv, M1n and array weights, serially and on 2 ranks. `"average"`
uses the mean over those directions: Gauss quadrature of the Derham for
array weights, 16-point Gauss–Legendre for callables.
- Array weights are reduced with a single `Allreduce`, instead of a
subcommunicator plus `Bcast`.
- The mass matrix is located in the composed operator by identity, and
the approximate inverse is applied exactly there. Both preconditioners
share this code (`_apply_composed`). Before, a mass matrix in the
left-most position was applied instead of inverted.
- `FFTSolver` copies its column, which was a view of the dense matrix
also used for the stencil factor. It now stabilises a singular matrix
once, at construction. `is_circulant` is vectorized.

PCG iterations (tolerance 1e-8) on HollowTorus, grid 12×16×6, degree
(2,3,2), as a reference for choosing `weight_reduction`:

| | M0 | M1 | M2 | M3 | Mv | M1n | M2n | Mvn |
|---|---|---|---|---|---|---|---|---|
| midpoint | 13 | 22 | 22 | 13 | 31 | 673 | 600 | 1000 (not converged) |
| average | 13 | 15 | 15 | 13 | 30 | 798 | 727 | 1000 (not converged) |

On HollowCylinder both need 2 iterations for every matrix.

### Kronecker preconditioners: base class, stiffness preconditioner

**`KroneckerPreconditioner`** is a new base class for preconditioners of
the form `P = B E · S Ã⁻¹ S · Eᵀ Bᵀ`:
- `Ã` (`matrix`) approximates the core operator, and its exact inverse
is `solver`. Subclasses build the approximation.
- `S` is an optional diagonal scaling `D̂^{1/2} D^{-1/2}`
(Loli–Sangalli–Tani), with `D` the diagonal of the core operator and
`D̂` that of the approximation.
- The base class holds the common interface (`solve`, `dot`, the
composition with `B E`, the properties).

**Mass preconditioners**
- `MassMatrixPreconditioner` is a subclass. It gets the options
`diagonal_scaling` (default `False`) and `dim_reduce=None` (unit weights
in all directions), and the method `update_mass_operator`. With
`diagonal_scaling=True`, the PCG iterations for M1 on a HollowTorus drop
from 20 to 8.
- `MassMatrixDiagonalPreconditioner` is now
`MassMatrixPreconditioner(dim_reduce=None, diagonal_scaling=True)`, kept
for its name in the solver options. It no longer assembles a separate 3d
logical mass matrix; `D̂` comes from the Kronecker approximation.
Results are unchanged up to rounding (1e-13), checked on HollowTorus and
IGAPolarCylinder, including `update_mass_operator`, serially and on 2
ranks.

**`StiffnessPreconditioner`** (new) preconditions the stabilized
stiffness operators `Gᵀ M1 G + σ M0`, `Cᵀ M2 C + σ M1` and `Dᵀ M3 D + σ
M2`:
- With Kronecker mass matrices, each diagonal block is a sum of
Kronecker products of 1d stiffness and mass matrices. It is inverted
exactly with `KroneckerSumSolver`, one solver per component (the
components have different spline types, e.g. DNN/NDN/NND for M1). For
grad, this is exact on the logical cube.
- `curl` and `div` neglect the off-diagonal blocks (block Jacobi). Block
Jacobi alone overestimates the operator on the kernel of the derivative,
so a kernel correction `σ⁻¹ d₋ P₋ d₋ᵀ` is added (default, needs `σ >
0`): with `d₋ = G` and the grad preconditioner for curl, and `d₋ = C`
and the curl block Jacobi for div.
- Weights: `weights="average"` (default) approximates the mass-matrix
weights by one common separable shape per direction plus a mean per
term, which keeps the fast diagonalization exact. `weights="unit"` uses
the logical cube. Diagonal scaling (default on) adds the geometry
pointwise; for it, the diagonal of the stiffness operator is computed by
probing.
- Polar splines are not supported yet.

PCG iterations (tolerance 1e-8), HollowTorus, grid 12×16×6, degree
(2,3,2):

| | no preconditioner | mass preconditioner | `StiffnessPreconditioner`
(defaults) |
|---|---|---|---|
| grad, σ = 0 | 187 | – | 41 |
| curl, σ = 1 | > 3000 | 716 | 46 |
| div, σ = 1 | > 3000 | 636 | 48 |

On a stretched cuboid (constant, anisotropic weights), grad needs 2
iterations. For large σ (1e4), the mass preconditioner is as good or
slightly better (about 13 vs 25 iterations).

## Tests

Run locally with the submodule branch: all preconditioner tests
(`test_preconditioner_transpose.py`, the new
`test_stiffness_preconditioner.py`,
`test_mass_matrices.py::test_mass_preconditioner*`,
`test_reduced_weight_1d`) pass serially and on 2 ranks (52 passed each).
The new stiffness tests cover the approximation against the operator on
the unit cube, PCG iteration counts for grad/curl/div, the kernel
correction, `weights="average"`, and independence of the MPI
decomposition. The rest is left to CI.

Don't raise the `feectools` pin in `pyproject.toml` until 0.7.0 is on
PyPI. Before merging, the submodule pointer should move to the merged
`devel-tiny` commit.

🤖 Generated with [Claude Code](https://claude.com/claude-code)

---------

Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Reviewers

Copilot code review Copilot
Copilot review effort
Applies to this pull request for everyone.Learn more about Copilot code review.
Copilot left review comments
@max-models max-models Awaiting requested review from max-models

Assignees

No one assigned

Labels

None yet

Projects

None yet

Milestone

No milestone

Development

Successfully merging this pull request may close these issues.

2 participants

@spossann