{"id":"410371f2-5914-4e53-b04a-c2a91e1fe1fa","arxiv_id":"2504.13821","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":5.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":0,"one_line_summary":"A few hundred lines of Julia deliver recursive TRMM and TRSM kernels that perform on par with cuBLAS and rocBLAS on large matrices and work on Apple Silicon for the first time.","lead":"This paper shows Julia code that runs triangular matrix multiply and solve on NVIDIA, AMD, and Apple GPUs from one portable API, reaching speeds close to vendor libraries on large matrices. It is the first report of these recursive kernels running on Apple Silicon.","discovery_kind":"new_application","skeptic_critique":{"model":"deepseek-v4-flash","headline":"GEMM backend is undisclosed, so the reported parity could be inherited from vendor BLAS rather than from the portable Julia layer; the portability claim is uninterpretable until this is identified.","rationale":"The paper's central claim is that a unified Julia implementation can be both portable and competitive. The recursive structure is conventional (KBLAS), and the paper explicitly credits that line of work, so the mathematical algorithm is not the issue. The load-bearing condition for the portability claim is that the expensive GEMM component is itself performance-portable, or at least that the authors disclose and scope the dependency on vendor GEMM. The paper contains no mention of which GEMM is used on any platform, no internal profiling to show GEMM dominates the runtime, and no error bars or repetitions in Fig. 3. The first omission is the most critical: it determines how to read every ratio in Fig. 3. This is a reproducibility and interpretability defect, not a correctness defect; the code is ostensibly available and the benchmarks could be rerun. The reader's conditional verdict is therefore appropriate, and the needed condition is precisely the missing GEMM-backend disclosure. No change to the verdict is needed beyond emphasizing that the revision must identify the GEMM backend, add repetitions/error bars, and ideally add a GEMM-only timing column.","tokens_in":6882,"tokens_out":7499,"duration_ms":72742,"concrete_test":"Inspect the code at the Zenodo DOI in ref. [3] and trace the GEMM dispatch for CuArray, ROCArray, and MTLArray (e.g., via `@which mul!(C, A, B)` or by grepping for CUBLAS/ROCBLAS/Metal calls). Then, on one NVIDIA and one AMD GPU, rerun the Fig. 3 benchmark in two configurations: (a) as published, and (b) with the GEMM call forced to a generic KernelAbstractions/GPUArrays GEMM. Also record the fraction of total runtime spent inside the GEMM call (e.g., with CUDA events wrapping each `mul!`). If configuration (b) collapses the ratio below parity, or if the GEMM fraction is near 100%, the headline result is inherited from vendor GEMM and should be restated as a demonstration of recursive orchestration, not of portable kernel performance.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The recursive scheme in §3.1 and §3.3 intentionally places almost all floating-point work in GEMM calls; only small base cases are handled by the custom kernels. The measured end-to-end times in Fig. 3 therefore reflect the performance of whatever `mul!`/GEMM routine is selected for each backend, and this is never stated. On NVIDIA/AMD the natural dispatch is CUBLAS/ROCBLAS GEMM; on Apple it may be a Metal Performance Shaders GEMM or a generic GPUArrays fallback. If the first two are true, the paper does not demonstrate a portable TRMM/TRSM at the kernel level; it demonstrates a portable recursive wrapper around vendor GEMMs, and the 'few hundred lines' claim omits the vendor libraries that supply most of the performance. If instead a generic Julia GEMM is used, that would be a much stronger result and needs to be stated explicitly. Without this information, the reader cannot tell whether the near-parity ratio in Fig. 3 is a property of the recursive algorithm, of the Julia abstraction, or of the underlying vendor GEMM.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","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].","tokens_in":7079,"tokens_out":6365,"duration_ms":62123,"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":[{"comment":"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.","section":"§3.1 and §3.3, Fig. 3"},{"comment":"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.","section":"§4.2 and §4.3"},{"comment":"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.","section":"§4.2"},{"comment":"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.","section":"§4"}],"minor_comments":[{"comment":"The caption says 'TRMM (bottom row)' where the surrounding text indicates the bottom row is TRSM; the caption should be corrected.","section":"Fig. 3 caption"},{"comment":"The contribution bullet mentions 'cuSOLVER and rocBLAS', but the experiments compare against cuBLAS and rocBLAS; the library names should be reconciled.","section":"§1 contribution bullet"},{"comment":"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.","section":"§3.3"},{"comment":"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.","section":"§4.2"},{"comment":"There are several typographical errors, including 'Silicion', 'recurive', and 'algorithmns'; the manuscript should be proofread before publication.","section":"Abstract and conclusion"}],"recommendation":"major_revision","confidential_remarks":"The core idea is sound and the code is public, but the performance claims cannot be evaluated without disclosing the GEMM backend and providing a full benchmarking protocol. The stress-test concern about undocumented GEMM dispatch is valid and load-bearing. The paper is otherwise within scope for a workshop on high-performance/heterogeneous computing and could be accepted after revision."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Two things to know. First, this is a real port, not a paper that pretends the recursion is new: it credits KBLAS for the recursive scheme and contributes a unified Julia implementation across CUDA, ROCm, and Metal, plus the first reported recursive TRMM/TRSM on Apple Silicon. Second, the headline “near-vendor performance” is plausible but under-supported as written, because the recursion delegates almost all floating-point work to GEMM and the paper never says which GEMM is called on each platform. Without that, the parity ratios in Fig. 3 could just be vendor BLAS behind a Julia facade.\n\nWhat the paper does well: the code is on Zenodo, the implementation is genuinely compact (a few hundred lines), and the authors are honest about inheriting the algorithm from KBLAS. Showing the same recursive structure running on three vendor GPUs with consistent trends is a legitimate engineering data point, especially for Apple Silicon where no vendor TRMM/TRSM exists to compare against. For a workshop paper, that is a solid contribution and the authors deserve credit for shipping reproducible code.\n\nThe soft spots are real but mostly fixable. The unknown GEMM backend is the one that matters: if the inner calls dispatch to cuBLAS and rocBLAS on those platforms, then the paper demonstrates a portable recursive wrapper around vendor GEMMs, not a portable kernel implementation. That is still useful, but it is a weaker claim than “unified Julia programs deliver near-vendor performance,” and the paper should say which one it is. The benchmarking protocol is also thin: no error bars, no raw timing tables, no mention of repetitions or warm-up, and the small-matrix TRSM gap is hand-waved as hardware idiosyncrasy. There is a minor inconsistency in the contributions list, which mentions cuSOLVER while the benchmarks use cuBLAS/rocBLAS. None of this is fatal, but the accepted version should pin down the GEMM dispatch and add basic benchmark metadata.\n\nBottom line: this is an honest, useful port that deserves a serious referee. The central engineering claim is plausible and the code is released; the missing details are the kind a workshop referee should request in one round, not a reason to desk-reject. I would send it to peer review, but I would not yet cite it as evidence that Julia itself delivers vendor-competitive triangular kernels until the GEMM backend is disclosed.","headline":"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.","tokens_in":7619,"tokens_out":1934,"would_cite":false,"duration_ms":20086,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":[],"pacs":[],"model":"deepseek-v4-flash","headline":"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…","keywords":["heterogeneous computing","task-based programming","recursive algorithms","Julia","KernelAbstractions","TRMM","TRSM","GPU performance portability"],"falsifier":"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.","tokens_in":6697,"feed_emoji":"⚡","tokens_out":9945,"duration_ms":86157,"temperature":0.7,"pith_summary":"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.","feed_headline":"One Julia codebase matches vendor GPU solvers on NVIDIA, AMD, Apple","feed_subtitle":"Recursive triangular kernels built on GEMM reach vendor-library speed at large sizes and debut on Apple Silicon","key_machinery":"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.","core_discovery":"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.","pith_inferences":["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."],"forward_implications":["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."],"supporting_citations":[{"why":"introduces the recursive TRMM/TRSM decomposition into GEMM calls that this implementation adapts.","marker":"[9]"},{"why":"extends the recursive triangular kernel framework across manycore architectures and supplies the base-kernel design.","marker":"[10]"},{"why":"is the released Julia implementation whose timing is benchmarked in the paper.","marker":"[3]"},{"why":"defines the BLAS Level 3 TRMM and TRSM operation signatures the paper implements.","marker":"[8]"},{"why":"provides the CUDA backend for NVIDIA GPU execution.","marker":"[7]"},{"why":"provides the AMDGPU backend for AMD GPU execution.","marker":"[17]"},{"why":"provides the Metal backend that makes Apple Silicon execution possible.","marker":"[6]"},{"why":"provides the oneAPI backend for the same kernel abstraction on Intel-family devices.","marker":"[5]"}],"fun_headline_variants":["Julia recursive kernels hit vendor speeds on all major GPUs","One Julia kernel matches cuBLAS on NVIDIA, AMD, Apple","Portable Julia code rivals vendor GPU libraries for triangular solves","Recursive Julia TRMM/TRSM matches vendor GPUs across all hardware","Julia's recursive triangular solvers near vendor speed on every GPU"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"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.","fun_headline_variants_meta":{"raw":{"variants":["Julia recursive kernels hit vendor speeds on all major GPUs","One Julia kernel matches cuBLAS on NVIDIA, AMD, Apple","Portable Julia code rivals vendor GPU libraries for triangular solves","Recursive Julia TRMM/TRSM matches vendor GPUs across all hardware","Julia's recursive triangular solvers near vendor speed on every GPU"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000486,"raw_usage":{"total_tokens":2369,"prompt_tokens":888,"completion_tokens":1481,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":504,"completion_tokens_details":{"reasoning_tokens":1394}},"tokens_in":504,"tokens_out":1481,"duration_ms":10583,"temperature":1.0,"reasoning_tokens":1394,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-16T11:58:38.962420+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"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.","supporting_citations":[{"cited_title":"A framework for dense triangular matrix kernels on various manycore architectures.Concurrency and Computation: Practice and Experience, 29 (22):e4187, 2017","cited_arxiv_id":null,"evidence_quote":"extends the recursive triangular kernel framework across manycore architectures and supplies the base-kernel design."},{"cited_title":"An Updated Set of Basic Linear Algebra Subprograms (BLAS)","cited_arxiv_id":null,"evidence_quote":"defines the BLAS Level 3 TRMM and TRSM operation signatures the paper implements."}],"review_version":1}