Pith. sign in

REVIEW 4 major objections 5 minor 4 cited by

Toward Portable GPU Performance: Julia Recursive Implementation of TRMM and TRSM

T0 review · 4 major / 5 minor · reviewed 2026-08-16 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗ Review JSON

Pith's one-line read A single recursive Julia implementation of TRMM and TRSM, built on GEMM calls and Julia's hardware-agnostic abstractions, runs on NVIDIA, AMD, and Apple Silicon GPUs and matches cuBLAS and rocBLAS performance for large matrices, while…

desk verdict A genuinely useful Julia port of KBLAS's recursive TRMM/TRSM with first Apple Silicon results, but the benchmark reporting hides the GEMM backend and noise, so the portability claim needs one more revision. read the letter →

arxiv 2504.13821 v1 pith:O3HCAAG5 submitted 2025-04-18 cs.MS cs.DC

Vicki Carrica , Maxwell Onyango , Rabab Alomairy , Evelyne Ringoot , James Schloss , Alan Edelman This is my paper · ORCID
classification cs.MScs.DC
keywords heterogeneouscomputingtask-basedprogrammingrecursivealgorithmsJuliaKernelAbstractionsTRMMTRSMGPUperformanceportability
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

The paper sets out to show that a unified Julia program can deliver near-vendor performance for triangular BLAS kernels across NVIDIA, AMD, and Apple GPUs. It restructures TRMM and TRSM recursively so that most floating-point work happens inside GEMM, the operation GPUs are best at, and uses Julia's multiple dispatch and metaprogramming, via the GPUArrays and KernelAbstractions packages, to keep one code path portable. For large matrices the reported runtimes are comparable to cuBLAS and rocBLAS, and the routines work on Apple Silicon, where no comparable vendor implementation existed. If the claim holds, performance portability would not require sacrificing speed for this central class of dense linear algebra kernels.

What carries the argument

The load-bearing mechanism is a recursive block decomposition that converts a triangular operation into mostly GEMM: for lower-triangular A split into diagonal block $A_{11}$, off-diagonal block $A_{21}$, and second diagonal block $A_{22}$, with B split into $B_1$ and $B_2$, the code first handles $A_{11}$ and $B_1$ recursively, then performs the update $B_2 = B_2 - A_{21} B_1$ (TRSM) or $B_1 = A_{21}^{\mathsf{T}} B_2 + B_1$ (TRMM) as a GEMM, then recurses on $A_{22}$ and $B_2$. This keeps the large compute-bound work in highly optimized GEMM kernels, isolates triangular dependencies in small base tiles, and improves memory reuse. Julia's multiple dispatch selects the correct kernel variant at each level based on side, triangularity, transpose, and solve-versus-multiply, so one code path covers all variants and backends.

What would settle it

Replace the GEMM calls inside the recursive kernels on each GPU with a deliberately slow but correct GEMM and rerun the same benchmarks; if the near-parity ratios survive, the recursion itself carries the performance, and if they collapse, the parity is inherited from the underlying GEMM rather than from the portable kernel.

Watch

Extended reading notes

Core claim

The central claim is that a hardware-agnostic recursive implementation of TRMM and TRSM in Julia reaches vendor-level throughput on large problems. The algorithm splits the triangular matrix into diagonal triangular blocks and off-diagonal blocks, recursively solves the top block, applies the off-diagonal update as a GEMM, and recursively solves the remaining block; below a tile threshold it invokes small base kernels. Because GEMM dominates, the kernels ride on the GPU's compute-bound strength and avoid triangular bottlenecks. Benchmarks report TRMM on rectangular inputs consistently faster than cuBLAS and rocBLAS, square TRMM within 50-200 percent of cuBLAS and at least 90 percent of rocBLAS, and TRSM at or near parity for larger sizes while matching at least two-thirds in square cases. The authors take this as evidence that unified, hardware-agnostic Julia abstractions can support production-grade level-3 BLAS kernels.

Load-bearing premise

The parity claim presumes the GEMM calls inside the recursion are already fast on each GPU and that the single-run timings are stable, but the paper does not identify which GEMM implementation is used or report error bars.

Editorial extensions

If this is right

  • A single Julia API can serve as a portable replacement for triangular BLAS on NVIDIA, AMD, and Apple GPUs, eliminating the need for vendor-specific kernels for these operations.
  • Large TRMM and TRSM workloads on Apple Silicon become possible for the first time, enabling triangular-solve-dependent algorithms on that hardware.
  • The recursive GEMM-centric structure extends naturally to upper-triangular and transposed variants, so new matrix layouts can be added by writing dispatch cases rather than new kernels.
  • The code's small size (a few hundred lines) lowers the maintenance burden of tracking three GPU vendor libraries.
  • At matrix sizes whose runtimes are below roughly ten milliseconds, the portable implementation can trail a vendor library, but the paper argues this regime is not where application time is dominated.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • A natural next experiment is to identify and benchmark the GEMM underlying the recursion on each platform; subtracting a pure GEMM baseline would reveal how much overhead the recursive wrapper actually adds.
  • The same block decomposition could be applied to other dense factorizations with triangular dependencies, such as Cholesky or LU updates, where GEMM-heavy rectangular updates and small triangular solves play a parallel role; the paper does not test this.
  • If Apple later releases an optimized vendor BLAS, this Julia implementation could serve as a portable baseline, and the measured gap on Apple would quantify how much performance is lost to hardware-agnostic abstractions rather than recovered via GEMM delegation.
  • The reported parity rests on single timing runs without variance; repeating the benchmarks with repetitions would show whether the remaining differences are real or within measurement noise.
Share X Bluesky LinkedIn Reddit HN

Signed reviews

No signed human review yet.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

4 major / 5 minor

Summary. The paper presents a Julia implementation of recursive TRMM and TRSM built on GPUArrays.jl and KernelAbstractions.jl. The recursive scheme decomposes triangular operations into GEMM calls interspersed with small base-kernel solves, and the paper claims that this hardware-agnostic implementation achieves performance comparable to cuBLAS and rocBLAS on large matrices while providing the first TRMM/TRSM support on Apple Silicon. Benchmarks are reported for NVIDIA A100, AMD MI100, and Apple M1 Pro platforms, with Fig. 3 showing runtime ratios against cuBLAS/rocBLAS. The authors state that the entire implementation is only a few hundred lines of code and make it publicly available [3].

Significance. If the performance parity claim holds, the paper is a useful case study in hardware portability: it shows that Julia's multiple dispatch, GPUArrays, and KernelAbstractions can present a single API for TRMM/TRSM across NVIDIA, AMD, and Apple GPUs with a small code footprint. The recursive algorithmic idea is not new, but the portability layer and the Apple Silicon demonstration are original. The code is public, the algorithmic description is clear, and this is a benchmark rather than a fitting paper, so circularity is not a concern. The main uncertainties are whether the measured parity is inherited from undisclosed vendor GEMM calls and whether the benchmarking protocol supports the claimed comparisons.

major comments (4)
  1. [§3.1 and §3.3, Fig. 3] The recursive scheme places almost all floating-point work in GEMM calls, but the paper never identifies which GEMM implementation is used on each platform. If `mul!` dispatches to cuBLAS/rocBLAS or to a Metal-specific GEMM, the end-to-end parity in Fig. 3 could be inherited from vendor GEMM, and the 'few hundred lines' claim would omit that dependency. The authors must state the exact GEMM routine per backend, and ideally include a GEMM-only baseline or a comparison against a generic Julia GEMM so that the recursive TRMM/TRSM overhead is isolated.
  2. [§4.2 and §4.3] The benchmark section reports only runtime ratios and provides no error bars, raw timing tables, repetition counts, warm-up procedures, or software-version details. Without these, the reader cannot distinguish the reported parity from run-to-run noise. The authors should supply raw timings, variance across repeated runs, and a precise description of the benchmarking protocol.
  3. [§4.2] The small-matrix TRSM underperformance is explained by 'hardware idiosyncracy' and 'algorithmic differences in the base kernel' without supporting data. Since this regime is clearly visible in Fig. 3, these statements are currently unsupported; the authors should either provide a kernel-level breakdown or explicitly label these comments as hypotheses.
  4. [§4] Section 4 reports only timings and performance ratios; no numerical correctness checks are presented for the Julia TRMM/TRSM implementations. For a numerical-library benchmark, performance comparisons are only meaningful if correctness is verified (e.g., residual norms for TRSM, error norms for TRMM). The authors should report accuracy checks or state that correctness is validated by existing test suites.
minor comments (5)
  1. [Fig. 3 caption] The caption says 'TRMM (bottom row)' where the surrounding text indicates the bottom row is TRSM; the caption should be corrected.
  2. [§1 contribution bullet] The contribution bullet mentions 'cuSOLVER and rocBLAS', but the experiments compare against cuBLAS and rocBLAS; the library names should be reconciled.
  3. [§3.3] The 'Memory optimization' bullet is too vague: 'shared memory and contiguous memory striding' does not describe the tile sizes, synchronization, or data layout used in the base kernels.
  4. [§4.2] The statement that Julia 'matches at least 2/3 of rocBLAS performance' is vague without a matrix-size threshold or a raw timing reference; please specify the sizes to which the statement applies.
  5. [Abstract and conclusion] There are several typographical errors, including 'Silicion', 'recurive', and 'algorithmns'; the manuscript should be proofread before publication.

Circularity Check

0 steps flagged · score 0.0 of 10

No circular derivation found; the benchmark parity claim is measured against external vendor libraries and does not reduce to fitted inputs or self-citations.

full rationale

This paper is an empirical porting and benchmarking study rather than a derivation or fitting exercise. The recursive TRMM/TRSM structure is attributed to prior KBLAS work by different authors, and the claimed contribution is a Julia implementation on NVIDIA, AMD, and Apple GPUs with measured runtime comparisons against cuBLAS and rocBLAS. No parameter is fitted to the benchmark data, no result is predicted from a model that embeds the same data, and no load-bearing argument is justified solely by a self-citation. The recursive decomposition in Section 3.1 is an algebraic block recursion whose correctness does not depend on the measured outcomes. The strongest adjacent concern is that the implementation delegates most floating-point work to GEMM calls and does not state which GEMM backend is selected on each platform; this is a transparency and attribution issue that affects how much of the observed parity comes from the portable Julia wrapper rather than an underlying vendor GEMM, but it is not circular reasoning because the GEMM performance is an external input, not an output of the paper's own derivation. The benchmark figures compare against external vendor libraries, so the central parity claim remains externally falsifiable. Therefore no circularity step is identified.

Assumptions & free parameters 0 free parameters · 3 assumptions · 0 invented entities

No free parameters or invented entities appear; this is an engineering benchmark paper. The main unstated inputs are the GEMM backend identity, package versions, and the benchmark protocol.

assumptions (3)
  • standard math The block recursive decomposition of triangular matrices into A11, A21, A22 and B1, B2 preserves the TRMM and TRSM identities.
    Section 3.1 partitions the matrices and applies the standard block triangular algebra; this is textbook linear algebra and the paper does not prove it, which is acceptable for this context.
  • domain assumption A high-performance GEMM routine is available and callable from Julia on NVIDIA, AMD, and Apple GPUs.
    The recursive scheme converts most work into GEMM calls (Section 3.1). The reported parity with cuBLAS/rocBLAS depends on those GEMM calls being fast on each platform, but the paper never names the GEMM backend.
  • domain assumption The benchmark timings are representative without repeated trials or error bars.
    Section 4 reports timing ratios in Figures 2 and 3 with no statistical treatment, implicitly assuming the numbers are stable and reproducible.

how reviews work

0 comments
Cite this review

Pith. "Pith review of Toward Portable GPU Performance: Julia Recursive Implementation of TRMM and TRSM." pith.science (2026). https://pith.science/paper/O3HCAAG5

@misc{pith2026250413821,
  author       = {Pith},
  title        = {Pith review of: Toward Portable GPU Performance: Julia Recursive Implementation of TRMM and TRSM},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/O3HCAAG5}},
  note         = {Machine review of arXiv:2504.13821}
}
read the original abstract

This paper presents a performant and portable recursive implementation of triangular matrix-matrix multiplication (TRMM) and triangular solve (TRSM) in Julia for GPUs, two kernels that underlie many linear-algebra algorithms. We restructure TRMM and TRSM so that most work is executed as general matrix-matrix multiplication (GEMM), improving use of the GPU memory hierarchy and reducing latency. Exploiting Julia's multiple dispatch and metaprogramming together with the GPUArrays and KernelAbstractions frameworks, we expose a single hardware-agnostic API that runs on NVIDIA, AMD, and Apple Silicon GPUs. For large matrices the recursive code reaches throughput comparable to vendor libraries such as cuBLAS and rocBLAS, while providing these routines on Apple Silicon for the first time. The entire implementation is only a few hundred lines of code, showing that unified Julia programs can deliver near-vendor performance across heterogeneous architectures.

Figures

Figures reproduced from arXiv: 2504.13821 by the authors.

Figure 1
Figure 1. TRMM/TRSM Recursive Illustration. In [PITH_FULL_IMAGE:figures/full_fig_p004_1.png] view at source ↗
Figure 2
Figure 2. Runtime of recursive unified TRMM (top row) and TRSM (bottom row) functions across different GPU hardware platforms (Apple, AMD, NVIDIA) as a function of the size of the matrix A ∈ R n×n for a rectangular matrix B ∈ R n×256, both of single precision. The figure shows a similar performance trend across hardware, demonstrating similar performance trends on three different hardware setups. performance differences could… view at source ↗
Figure 3
Figure 3. Runtime ratio of cuBLAS/rocBLAS versus the Julia implementation of TRMM (top row) and TRMM (bottom row) in function of the size of matrix A ∈ R n×n . Higher values indicate that the Julia implementation is faster, 100% indicates equal performance. The left two figures are for a matrix ∈ R n×256 having a set width. The right two figures show the case of a square matrix B ∈ R n×n . The figures demonstrates the unified… view at source ↗

Discussion (0). Continue with ORCID to comment.

Forward citations

Cited by 4 Pith papers

Reviewed papers in the Pith corpus that reference this work. Sorted by Pith novelty score. Full citation record

  1. Hierarchical Recursive Precision for Accelerating Symmetric Linear Solves on MXUs

    cs.DC 2026-01 conditional novelty 6.0 of 10

    A tree-structured recursive Cholesky solver assigns FP16 to off-diagonal blocks and higher precision to diagonal blocks, achieving over 5x speedup on NVIDIA H200 and AMD MI300X with better accuracy than pure half precision.

  2. Accelerating Bidiagonalization of Banded Matrices through Memory-Aware Bulge-Chasing on GPUs

    cs.DC 2025-10 conditional novelty 6.0 of 10

    A memory-aware GPU bulge-chasing algorithm reduces banded matrices to bidiagonal form, achieving >100x speedups over CPU libraries at 32k sizes.

  3. Performant Unified GPU Kernels for Portable Singular Value Computation Across Hardware and Precision

    cs.DC 2025-08 conditional novelty 6.0 of 10

    A unified Julia implementation of two-stage QR SVD achieves near-cuSOLVER performance across four GPU vendors and three precisions, including firsts for Apple Metal and half precision.

  4. Julia for CFD: A Critical Survey of Ecosystem, Performance, and Composability

    cs.CE 2026-08 conditional novelty 3.0 of 10

    Julia is now credible for several research-facing CFD regimes, but its advantage is integration and composability, not unique performance, and it is not yet a complete industrial CFD platform.

Reference graph

Works this paper leans on

17 extracted references · 10 canonical work pages · cited by 4 Pith papers

  1. [3]

    NextLA.jl: Next-Gen Linear Algebra

    Rabab Alomairy, Evelyne Ringoot, Sophie Xuan, Vicki Carrica, Maxwell Onyango, and Julian Samaroo. NextLA.jl: Next-Gen Linear Algebra. Zenodo, https://doi.org/10.5281/zenodo.15049222, 2025

    doi:10.5281/zenodo.15049222 2025

  2. [1]

    Communication Avoiding {LU} with Tournament Pivoting in{SLATE},{SWAN} No

    Rabab Alomairy, Mark Gates, Sebastien Cayrols, and Dalal Sukkari. Communication Avoiding {LU} with Tournament Pivoting in{SLATE},{SWAN} No. 18. 2022

    work page 2022

  3. [2]

    DynamicTaskSchedulingwith Data Dependency Awareness Using Julia

    RababAlomairy,FelipeTome,JulianSamaroo,andAlanEdelman. DynamicTaskSchedulingwith Data Dependency Awareness Using Julia. In2024 IEEE High Performance Extreme Computing Conference (HPEC). IEEE, 2024

    work page 2024

  4. [4]

    High-Performance Scientific Applications Using Mixed Precision and Low- Rank Approximation Powered by Task-based Runtime Systems

    Rabab M Alomairy. High-Performance Scientific Applications Using Mixed Precision and Low- Rank Approximation Powered by Task-based Runtime Systems. 2022

    work page 2022

  5. [5]

    oneAPI.jl, January 2025

    Tim Besard. oneAPI.jl, January 2025. URLhttps://doi.org/10.5281/zenodo.14615352

    doi:10.5281/zenodo.14615352 2025

  6. [7]

    Effective Extensible Programming: Unleash- ing Julia on GPUs.IEEE Transactions on Parallel and Distributed Systems, 30(4):827–841, 2019

    Tim Besard, Christophe Foket, and Bjorn De Sutter. Effective Extensible Programming: Unleash- ing Julia on GPUs.IEEE Transactions on Parallel and Distributed Systems, 30(4):827–841, 2019. https://doi.org/10.1109/TPDS.2018.2872064

    arXiv 2019

  7. [8]

    An Updated Set of Basic Linear Algebra Subprograms (BLAS)

    L Susan Blackford, Antoine Petitet, Roldan Pozo, Karin Remington, R Clint Whaley, James Demmel, Jack Dongarra, Iain Duff, Sven Hammarling, Greg Henry, et al. An Updated Set of Basic Linear Algebra Subprograms (BLAS). ACM Transactions on Mathematical Software, 28 (2):135–151, 2002

    work page 2002

  8. [9]

    Redesigning Triangular Dense Matrix Com- putations on GPUs

    Ali Charara, Hatem Ltaief, and David E Keyes. Redesigning Triangular Dense Matrix Com- putations on GPUs. In Euro-Par 2016: Parallel Processing, pages 477–489. Springer, 2016. https://doi.org/10.1007/978-3-319-43659-3_35

    doi:10.1007/978-3-319-43659-3_35 2016

Show all 17 references
  1. [10]

    A framework for dense triangular matrix kernels on various manycore architectures.Concurrency and Computation: Practice and Experience, 29 (22):e4187, 2017

    Ali Charara, David Keyes, and Hatem Ltaief. A framework for dense triangular matrix kernels on various manycore architectures.Concurrency and Computation: Practice and Experience, 29 (22):e4187, 2017. https://doi.org/10.1002/cpe.4187

    2017 doi

  2. [11]

    Valentin Churavy.Language Evolution for Parallel and Scientific Computing. Ph.d. thesis, Mas- sachusetts Institute of Technology, Department of Electrical Engineering and Computer Science, Cambridge, MA, September 2024. Licensed under a CC BY-NC-ND 4.0 license

    2024

  3. [12]

    Bridging HPC Communities through the Julia Programming Language

    Valentin Churavy, William F Godoy, Carsten Bauer, Hendrik Ranocha, Michael Schlottke- Lakemper, Ludovic Räss, Johannes Blaschke, Mosè Giordano, Erik Schnetter, Samuel Omlin, et al. Bridging HPC Communities through the Julia Programming Language. arXiv preprint arXiv:2211.02740, 2022

    2022 arXiv

  4. [13]

    An Evaluative Comparison of Performance Portability across GPU Programming Models.arXiv preprint arXiv:2402.08950, 2024

    Joshua H Davis, Pranav Sivaraman, Isaac Minn, Konstantinos Parasyris, Harshitha Menon, Gior- gis Georgakoudis, and Abhinav Bhatele. An Evaluative Comparison of Performance Portability across GPU Programming Models.arXiv preprint arXiv:2402.08950, 2024

    2024 arXiv

  5. [14]

    Programming Heterogeneous Architectures Uing Hierarchical Tasks.Concurrency and Computation: Practice and Experience, 35(25):e7811, 2023

    Mathieu Faverge, Nathalie Furmento, Abdou Guermouche, Gwenolé Lucas, Raymond Namyst, Samuel Thibault, and Pierre-andré Wacrenier. Programming Heterogeneous Architectures Uing Hierarchical Tasks.Concurrency and Computation: Practice and Experience, 35(25):e7811, 2023

    2023

  6. [15]

    Evo- lution of the SLATE Linear Algebra Library.The International Journal of High Performance Computing Applications, 39(1):3–17, 2025

    Mark Gates, Ahmad Abdelfattah, Kadir Akbudak, Mohammed Al Farhan, Rabab Alomairy, Daniel Bielich, Treece Burgess, Sébastien Cayrols, Neil Lindquist, Dalal Sukkari, et al. Evo- lution of the SLATE Linear Algebra Library.The International Journal of High Performance Computing Ap...

    2025

  7. [16]

    Dynamic Task Discovery in PaRSEC: A Data-Flow Task-Based Runtime

    Reazul Hoque, Thomas Herault, George Bosilca, and Jack Dongarra. Dynamic Task Discovery in PaRSEC: A Data-Flow Task-Based Runtime. InProceedings of the 8th Workshop on Latest Advances in Scalable Algorithms for Large-Scale Systems, pages 1–8, 2017

    2017

  8. [17]

    Julian Samaroo, Anton Smirnov, Valentin Churavy, Ludovic Räss, Torrance Hodgson, Alexis Mon- toison, Wiktor Phillips, Ali Ramadhan, Jason Barmparesos, Tim Besard, Julia TagBot, Michel 10 Carrica et al. Schanen, Carsten Bauer, Mosè Giordano, Takafumi Arakaki, Stephan Antholzer,...

    2023 doi

  9. [18]

    Synthesizing Numerical Linear Algebra using Julia

    Sophie Xuan, Evelyne Ringoot, Rabab Alomairy, Felipe Tome, Julian Samaroo, and Alan Edel- man. Synthesizing Numerical Linear Algebra using Julia. In 2024 IEEE High Performance Extreme Computing Conference (HPEC). IEEE, 2024

    2024

Pith tools

Reviewed August 16, 2026 · model on record in the stance chip above.