Add experimental sparse matrix-matrix multiplication

[liaoweiyang2017]

Summary

This draft adds experimental spmm support following the SpMV layout: one public module, five format-specific submodules, and a shared private implementation submodule. The array-kernel revision is in 53a9bb70; 3d1e5158 only relocates the changelog entry to avoid an upstream conflict.

• Sparse/dense and dense/sparse forms retain the existing SpMV scaling and operation semantics.
• Sparse/sparse forms accept same-format, same-kind COO, CSR, CSC, ELL and SELLC matrices with sparse_full storage.
• Public spmm_prepare builds C independently of numerical values, with independent N/T/H operations on A and B.
• Five public spmm_kernel_* interfaces operate on plain arrays and an integer workspace. No persistent plan or structural snapshots are used.
• High-level spmm calls preparation by default; allow_resize=.false. reuses the supplied result structure.

type(CSR_dp_type) :: a, b, c
integer(ilp), allocatable :: work(:)

! Initialize a and b first.
call spmm_prepare(a,b,c)
allocate(work(c%ncols))
call spmm(a,b,c,allow_resize=.false.,work=work)
! Change only input values, then reuse the same result and workspace.
call spmm(a,b,c,allow_resize=.false.,work=work)
! Rebuild after a structural change when necessary.
call spmm(a,b,c)

The caller must ensure that C contains the complete product pattern; a valid superset is accepted. Structural changes are not detected automatically. The high-level wrapper checks operations, dimensions and storage/workspace extents; raw kernels assume conforming arrays. The result must not alias an input, and concurrent calls require separate result/work buffers. Stored zeros and cancellation do not remove prepared positions.

Local validation

• 17/17 sparse tests pass with GNU 13.4 runtime checking (real/complex SP, DP and QP), GNU 16 runtime checking (SP/DP), and GNU 13 FPM Release. Tests cover all five formats and nine N/T/H pairs, rectangular/empty products, scaling, unordered COO inputs, valid result supersets, storage reuse and input errors.
• Both examples compile and run; FORD generates the new interfaces and specification.
• Allocation probes record zero allocations for supplied-work N×N CSR/CSC/ELL/SELLC kernels and wrappers. COO grouping and transposed paths still allocate temporary integer views; C is retained.
• Six large fixtures (10,000/30,000/50,000 square, 10/30 entries per column) yielded 36 full pointer/index/value comparisons against the previous implementation and Julia. No mismatches at rtol=atol=1e-12; maximum absolute difference 4.44e-16.

Performance evaluation

Apple M5 Pro, GNU Fortran 13.4.0 -O3, Float64/Int32 CSC, 5 warm-ups, 15 samples and three independent processes. At 50,000 square with 30 entries per input column:

Stage	Previous plan design	Array design
Preparation	95.7 ms	71.7 ms
Repeated numerical kernel	61.6 ms	49.9 ms
Complete call rebuilding C	159.6 ms	115.7 ms

For the same fixture, the original fused Fortran implementation takes 92.9 ms for a complete call. Julia's official allocating multiplication takes 160.0 ms; a separate local Julia array reference takes 44.6 ms for numerical reuse (not an official SparseArrays API). At 50,000 square and 10 entries per column, the new Fortran complete call takes 13.2 ms versus Julia's 11.2 ms. These measurements cover one structured family on one machine, not a general ranking. Benchmark and allocator probes remain outside the PR diff.

The two-phase path improves repeated use but retains extra one-shot symbolic traversal and initialization. Whether to maintain a separate fused one-shot path remains open for discussion, as do the caller-owned pattern contract and same-format scope. Hosted CI for this revision requires maintainer approval; local passes do not establish hosted CI success.

AI assistance

OpenAI Codex assisted with design, implementation, tests, benchmark harnesses and documentation. The validation and timing results above come from actual local runs.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 68.12%. Comparing base (90f135a) to head (017a960).

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1215      +/-   ##
==========================================
- Coverage   68.20%   68.12%   -0.09%     
==========================================
  Files          19       19              
  Lines        2378     2378              
==========================================
- Hits         1622     1620       -2     
- Misses        756      758       +2     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[jalvesz]

Thanks for the proposal @liaoweiyang2017

Just a few first comments before a full review:

• I would propose making the construction of the target sparse object a public procedure that can used independently of the spmm kernel
• Following the previous comment, I would propose to allowing the use of the spmm kernel given a target matrix that is already pre-computed to proposer dimensions

Why? For performance sake, it might be desired to call the spmm kernel in an itereative manner without going through runtime reallocations of the target matrix but just to redo the operations with the same triplet of matrices.

In other scenarios, like computing the exponent of a matrix, indeed an incremental increase of the matrix sparsity pattern might be desired, so we need a way to control whether to just let target matrix to be resized or just recycle the same existing matrix.

I see you managed the combinations of different sparse types, I would have kept the PR reduced to same-type products, open for debate.

[liaoweiyang2017]

Thanks @jalvesz for the suggestions. I have pushed an update in commit 859096f9 to separate symbolic preparation from numerical multiplication:

• spmm_prepare(a,b,c,plan) constructs the target pattern and reusable integer workspace independently of numerical values.
• spmm_kernel(a,b,c,plan) updates only numerical values. It checks the input and target structures before modifying the target and does not allocate result/work buffers.
• spmm(a,b,c,plan=plan) reuses a compatible plan by default. allow_resize=.true. explicitly permits rebuilding when the structure changes. The call without a plan retains the convenience of preparing and computing a product.
• Sparse-by-sparse overloads now require the same storage format and numeric kind for both factors and the target. Mixed-format products have been removed from this draft. COO, CSR, CSC, ELL and SELLC remain supported.

The symbolic pattern includes stored zeros and positions whose numerical contributions cancel, so changes in numerical values do not force preparation again. A change to indices or storage layout requires a new preparation; dimensions alone are not sufficient to establish compatibility.

I added tests for numerical reuse, unchanged storage addresses, structure mismatches, explicit rebuilding, zero scaling, cancellation, invalid storage and empty products, and updated the examples and specification. The local sparse suite passes in FPM Debug and CMake Release, including the configured real/complex kinds. An allocator-counting probe also records zero allocations for repeated prepared calls in each of the five formats. A complete 50,000-by-50,000 CSC result agrees with Julia in both pattern and values.

The reusable path reduces repeated-call time on the local fixture. Preparation adds time and integer memory compared with the previous fused one-shot implementation; I will keep that distinction explicit in the validation notes. A separate comparison now covers six large CSC fixtures. The Fortran prepared path is close to a local Julia reference implementing the same checked numerical-reuse workflow, with about 1.4%–13.0% higher call times in these measurements. That Julia reference is not an official SparseArrays API. Comparisons with the official allocating A*A call are reported separately, and the one-shot performance depends on the sparsity pattern.

The names and plan-based interface are still experimental, and I would appreciate feedback on this separation before considering the PR ready for a full review.

[jalvesz]

Hi @liaoweiyang2017

A first general comment: for the sake of transparency, could you please mention which Gen AI system are you using to help you with the development?

Coming back to the subject. I have the impression that the current design uses OOP quite heavily and I wonder if it is truly needed? Some aspects do require but I believe that OOP should be kept at minimum, specially in performance sensitive areas of code where plain array manipulations could be enough.

I'm concerned about why the split of matrix creation and application of the actual kernel introduced a measurable loss in performance. I would have imagined:

• A non-oop kernel which just applies the spmm kernel using plain array syntax for the 3 matrices. This kernel assumes the 3 matrices to conform, a minimal check of high level dimensions could suffice
• A spmm_prepare procedure to prepare allocate and build C as a function of A, B and the operator on each (transpose or not).
• the high-level spmm with an additional boolean argument allow_resize (or something of the sorts), which would then call or not the prepare procedure and then the kernel.

Before modifying your current proposal, could you evaluate this comment and give a feedback on what it would imply in terms of both performance and code quality? Thanks!

[liaoweiyang2017]

Thanks @jalvesz for the feedback. I am using OpenAI Codex to assist with the design, implementation, tests and benchmark harnesses. The reported results come from actual local compiler/test runs, and I have added this disclosure to the PR description.

I evaluated the array-based design locally and have now pushed the tested candidate in 53a9bb70 so the implementation and trade-offs are inspectable. This remains a draft, and the interface is open for discussion. 3d1e5158 only relocates the changelog entry to avoid an upstream conflict.

I agree that the persistent plan was more machinery than this interface needs. The previous numerical path was not doing virtual dispatch for every multiplication; the main removable work was structural snapshots, full index comparisons on every call, identity slot maps for CSR/CSC, and repeated copying during preparation.

The revision has no custom plan/structure types or polymorphic dispatch in the SpMM implementation:

• Five spmm_kernel_* interfaces operate directly on the three matrices' arrays and an integer workspace, following the SpMV layout.
• spmm_prepare(A,B,C,op_a=...,op_b=...) constructs the structural result independently of numerical values. Each operand supports N/T/H.
• spmm(A,B,C,allow_resize=.false.,work=...) skips preparation and reuses C. The default true value calls preparation before the array kernel.

The contract is intentionally lighter: the high-level wrapper checks operations, dimensions, storage and buffer/workspace extents. The raw kernel assumes conforming storage and a C pattern containing every required output position. A valid precomputed superset is accepted. Changes to values alone need no preparation; changes to indices or operations require the caller to establish that C is still sufficient or call preparation again. Exact structural-change detection is no longer provided automatically.

On the same Apple M5 Pro with GNU Fortran 13.4.0, -O3, Float64/Int32 CSC inputs, 5 warm-ups, 15 samples and three independent processes, the 50,000-square, 30-entries-per-column case gave:

Stage	Previous plan design	Array revision
Preparation	95.7 ms	71.7 ms
Repeated numeric kernel	61.6 ms	49.9 ms
Complete call rebuilding C	159.6 ms	115.7 ms

The original fused implementation, rebuilt with the same compiler/options, took 92.9 ms for the complete call. Thus removing the extra state recovers substantial performance, but a separate symbolic traversal and initialization still have a measurable one-shot cost. Iterative use amortizes preparation. Preserving the best fused one-shot time as well would require an additional execution path and its maintenance cost.

For the same large case, Julia's official allocating multiplication took 160.0 ms. A separately labelled local Julia array reference took 44.6 ms for the numerical kernel; this is not an official SparseArrays API. At 50,000 square and 10 entries per column, the Fortran complete call took 13.2 ms versus Julia's 11.2 ms. These are results for one structured family on one machine, not a general speed ranking.

The direct N×N CSR/CSC/ELL/SELLC paths allocate nothing when workspace is supplied. COO grouping and transposed inputs still use transient integer views; they do not reallocate C or copy input numerical values. A strict allocation-free requirement for every format/operation would need an explicit larger workspace layout or stronger input-ordering requirements. For the large CSC example, the old plan's integer storage was about 302 MiB by array-size accounting; the direct numeric workspace is about 0.19 MiB, excluding the matrices and preparation temporaries.

Locally, the 17 sparse tests pass with GNU 13/16 runtime checking and GNU 13 FPM Release (GNU 13 includes real/complex SP/DP/QP; the GNU 16 build covers SP/DP). They cover all five formats, all N/T/H pairs, rectangular and empty cases, and result storage reuse. Both examples and the FORD build pass. Six large cases produced 36 complete pointer/index/value comparisons against the previous implementation and Julia, with no mismatches at rtol=atol=1e-12 (maximum absolute difference 4.44e-16). Hosted workflows currently await maintainer approval.

In code-quality terms, the persistent state and full-structure validation disappear, while responsibility for a reusable result pattern becomes explicit. The source is slightly longer because it adds independent operations and five public array interfaces; the high-level usage becomes simpler. I would appreciate your view on this contract and on whether the initial implementation should prioritize this two-phase reuse path or also retain a fused one-shot specialization.

[jalvesz]

Thanks for the detailed explanation and figures.

I guess that, given the numbers, we could reconsider coming back to a fused call and keep the boolean flag for avoiding the reallocations if needed. That way the C matrix would be formed once already with the first application and subsequent calls can be more economical by passing the boolean to .false.

What do you think about it @liaoweiyang2017 ?

Maybe @jvdp1 has a different take on this?