Stencil 3D device kernels: 6D matrix views instead of the 2p+1 assumption

[max-models]

Stack: part 2 of 3 — based on #96 (device-kronecker-solve, merge that first), next: #98 (inner-on-device).

device-kronecker-solve is merged into this branch (merge commit d16a6ad, no rebase, so review comments stay); the only conflict was CUDA_STRATEGY.md, where both new sections are kept ("Kronecker solver on the device (#96)" first, then "6D matrix views"). The diff against the base shows only this PR's change. When #96 is merged, GitHub retargets this PR to devel-tiny.

Companion struphy PR: struphy-hub/struphy#706 (part 2 of the struphy stack struphy-hub/struphy#703 → struphy-hub/struphy#706 → struphy-hub/struphy#716); its feectools submodule points at this branch's new head d16a6ad, which contains #96.

---

Corresponding PR in struphy: struphy-hub/struphy#706

Summary

The 3D device stencil kernels from #88 (stencil_dot_3d, stencil_transpose_3d) take the matrix data as 6D views (Array6D<double>, cunumpy >= 0.6.1) instead of raw pointers, and all dot/transpose kernels read the number of diagonals from the matrix data instead of assuming 2 * p + 1. StencilMatrix.dot, vdot and transpose therefore no longer raise NotImplementedError for matrices with fewer diagonals than the pads of their spaces allow (blocks between spaces of different degree, derivative-type stencils). Related to struphy-hub/struphy#650 and struphy-hub/struphy#688 (6D array views in cunumpy).

Changes

• CUDA, 3D: stencil_dot_3d(Array6D<double> mat, Array3D<double> x, Array3D<double> out, s_in, p_in, add, s_out, e_out, p_out) and stencil_transpose_3d(Array6D<double> mat, Array6D<double> matT, s_in, p_in, add, s_out, e_out, p_out). Strided views as in 1D/2D (Array2D/Array4D), so the hand-computed strides are gone. The pyccel versions take float[:, :, :, :, :, :] with the same arguments in the same order.
• Diagonals from the data (all six kernels, pyccel and CUDA, 1D to 3D): along each direction the matrix has n = mat.shape[ndim + k] diagonals, its pads are q = (n - 1) // 2, diagonal d of row i is the column i - q + d (at x[i - q + d - s_in + p_in]), and the last owned row uses n - 1 + add. The transpose maps diagonal d of matT to diagonal q + i - j of mat. With q = p this is the old loop, in the same order. The pyccel kernels are one loop nest that picks the number of diagonals per row, instead of the 2/4/8 spelled-out (interior, last row) combinations.
• e_in removed from stencil_transpose_1d/2d/3d. It was only needed for the row extents of the raw pointer in 3D, and the argument list is the same for all dimensions.
• StencilMatrix: the 2 * p + 1 check is gone from dot, vdot and transpose. transpose(out=...) asserts that out has the pads of the matrix. Spaces with shifts > 1 still raise NotImplementedError (see below).
• Tests:
  • MATRIX_CASES (GPU parity and CPU emulation) gets 2 new 1D cases, 2 in 2D and 4 in 3D. They cover fewer diagonals in some directions, non-periodic rectangular blocks between spaces of different size per direction, a derivative-type rectangular block with fewer diagonals in every direction, and pads 0 (one diagonal).
  • test_device_matvec.py: the old "reject" test is replaced by test_fewer_diagonals_match_dense_reference, which checks dot, vdot, transpose and transpose(out=...) against toarray() for 8 such matrices. A new test checks that shifts > 1 raise.
  • test_mpi_device.py: a matrix with fewer diagonals against the global field, and the adjoint identity of its transpose for a non-symmetric (one-sided) stencil.
• CUDA_STRATEGY.md: new section "6D matrix views". The 2 * p + 1 / raw-pointer limitation and the open question about 6D views are removed, and shifts > 1 are listed as an open question.
• Dependency: cunumpy >= 0.6.1 was already required on devel-tiny (Update to cunumpy 0.6.1 #94), so it is unchanged.

Behaviour changes

• Matrices with fewer diagonals than 2 * p + 1 (StencilMatrix(V, W, pads=q) with q < p) work in dot, vdot and transpose on both backends, where they used to raise.
• Matrices that already worked give bitwise the same results. I checked this by running the old and new pyccel kernels uncompiled on every old case.
• stencil_transpose_<n>d takes one argument fewer (e_in). Only StencilMatrix calls these kernels; struphy does not call them directly.
• Shifts > 1 still raise. The stencil kernels never handled them, and there is no reference to match. psydac's general kernels (matvec_<n>d, transpose_<n>d) disagree with toarray() for shifts > 1. Their transpose is also not the adjoint of their product, and psydac never tested them ("TODO: verify for s>1"). Supporting shifts would first need a verified data layout.
• NumPy timings are unchanged (32×32×16, degree 3, compiled C): dot 5.9 → 5.1 ms, transpose 6.8 → 6.8 ms.

Testing

Local runs on macOS without a GPU (cunumpy 0.6.1, pyccel 2.2.1 with C, Open MPI 5). The GPU tests were not run. Before → after:

• test_cuda_parity.py, test_cuda_emulation.py, test_device_matvec.py: 29 passed / 31 skipped → 37 passed / 49 skipped. More GPU parity cases are skipped now, and the CPU emulation covers all cases, the new ones included.
• Serial feectools/linalg -m "not mpi and not petsc": 7942 passed → 7950 passed.
• Serial feectools -m "not mpi and not petsc": 9456 passed, 6 failed → 9464 passed, 6 failed. The 6 failures are the pre-existing ones in ddm/tests/test_cart_2d.py/test_cart_3d.py.
• mpirun -n 2 … feectools/linalg -m "mpi and not petsc" --with-mpi: 738 passed → 739 passed. test_mpi_device.py also passes with 1, 3 and 4 ranks.
• Extra checks (not committed):
  • CPU emulation of every matrix case compiled with -DCUNUMPY_BOUNDS_CHECK: no out-of-bounds index.
  • dot, vdot and transpose against toarray() for every case with a dense meaning.

🤖 Generated with Claude Code