linear algebra

This work addresses symmetric positive-definite block-tridiagonal linear systems arising in time-dependent estimation and optimal control. We propose a GPU-accelerated parallel solver based on recursive Schur complement reduction, which hierarchically decomposes the original problem into independent, batch-processable subproblems. To maximize GPU throughput, we design customized batched BLAS/LAPACK kernels and optimize task partitioning and memory access patterns via CUDA. The resulting open-source, cross-platform solver—TBD-GPU—demonstrates substantial speedups over state-of-the-art CPU-based sparse solvers (e.g., CHOLMOD, HSL MA57) across multiple benchmark problems, while matching the performance of NVIDIA cuDSS. These results validate the effectiveness and advancement of structure-aware batched computation for solving structured dense systems on GPUs.

GPU-accelerated triangular matrix–matrix multiplication (TRMM) and triangular solve (TRSM) kernels suffer from poor portability and inconsistent performance across heterogeneous GPU architectures (NVIDIA, AMD, Apple Silicon). Method: This paper proposes a unified algorithmic framework based on recursive divide-and-conquer and GEMM redirection. It reformulates TRMM/TRSM as sequences of standard GEMM calls, optimizes memory access patterns and compute–memory overlap, and leverages Julia’s multiple dispatch and metaprogramming to achieve hardware-agnostic abstraction. Contribution/Results: The implementation—requiring only ~500 lines of code—enables single-API deployment across all target platforms. It delivers the first high-performance TRMM/TRSM implementation on Apple Silicon GPUs; for large matrices, throughput approaches that of cuBLAS and rocBLAS, while cross-architecture performance variance is substantially reduced. This work breaks the long-standing trade-off between high performance and high portability in GPU linear algebra kernels.

This work addresses the lack of efficient solutions for computing a large batch of small-scale singular value decompositions (SVDs) on GPUs. The authors propose a GPU-accelerated batched SVD solver based on the one-sided Jacobi algorithm, co-designed with hardware architecture to exploit fine-grained parallelism, optimize memory access patterns, and support multiple floating-point precisions. Implemented on both NVIDIA and AMD GPU platforms, the solver demonstrates exceptional robustness and scalability across diverse matrix shapes, conditioning numbers, and precision configurations. Experimental results show that the proposed method significantly outperforms existing vendor-provided libraries and open-source solvers in terms of computational performance while maintaining numerical reliability.

This study addresses the performance bottleneck of large-scale square matrix multiplication (N ≥ 10,000) in deep learning and scientific computing. We systematically benchmark five mainstream parallel implementations—CuBLAS, native CUDA, CPU-based BLAS, OpenMP, and C++ standard threads—using a unified benchmarking framework, strictly controlled variables, and repeated trials under identical hardware and data conditions. Double-precision FLOPS serves as the primary empirical metric. Results demonstrate that CuBLAS significantly outperforms all four alternatives, with performance ranking exhibiting high statistical significance (p < 5×10⁻¹²). This work quantifies the practical efficiency disparities among parallel programming paradigms for ultra-large matrix multiplication and provides reproducible, evidence-based guidance for selecting high-performance linear algebra libraries. It further underscores the irreplaceable optimization advantages of hardware-aware libraries—particularly CuBLAS—in GPU-accelerated computation.

This work addresses the inefficiencies of conventional GPU linear algebra libraries, which often fail to fully exploit hardware capabilities due to redundant memory accesses and runtime overhead. To overcome these limitations, the authors propose Bandicoot, a toolkit that leverages C++ template metaprogramming to perform expression fusion at compile time, thereby automatically generating highly optimized GPU kernels that saturate memory bandwidth—without relying on just-in-time compilation or runtime scheduling. Bandicoot provides an API compatible with Armadillo, facilitating straightforward migration of existing CPU codebases. Experimental results demonstrate that Bandicoot significantly outperforms PyTorch, TensorFlow, and JAX across multiple benchmarks, achieving substantial speedups in several scenarios.

Existing LLM-driven approaches for GPU kernel optimization are largely confined to machine learning scenarios, lacking cross-domain generality and a systematic evaluation benchmark. To address this gap, this work proposes CUDAMaster—a hardware-aware, multi-agent automated optimization framework that integrates performance profiling with an automated compilation toolchain to generate CUDA kernels across diverse workloads under FP32 and BF16 precision. Furthermore, we introduce MSKernelBench, the first comprehensive benchmark encompassing algebraic operations, LLM operators, sparse matrix computations, and scientific computing. Experimental results demonstrate that CUDAMaster consistently outperforms Astra by an average of approximately 35% across most operators, with certain kernels matching or even surpassing the performance of the proprietary cuBLAS library.

This work presents the first systematic evaluation of Rust’s suitability for high-performance sparse linear algebra, addressing the longstanding trade-off between performance and memory safety in traditional scientific computing that relies on C/C++ and Fortran. The authors natively implement core operations—including sparse matrix-vector multiplication, the Lanczos Krylov method, and matrix exponential computation—leveraging compile-time monomorphization, SIMD vectorization, and careful FFI boundary analysis. Comprehensive benchmarks against established libraries such as Intel oneMKL, Eigen, PETSc, and PSBLAS demonstrate that Rust achieves performance on par with Eigen and PSBLAS in CSC format, approaching state-of-the-art levels while preserving memory safety. However, it still lags behind PETSc in block CSR optimizations, highlighting both the promise and current limitations of Rust for building efficient, safe numerical software stacks.

This work addresses the limitations of conventional sparse matrix formats on GPUs for sparse matrix-vector multiplication (SpMV), which suffer from high memory overhead, frequent data movement, and insufficient support for flexible precision. The authors propose a precision-agnostic sparse storage format built upon the SELL structure, integrating column index delta encoding with a single-word packing scheme that compresses both index deltas and numerical values into a single word. This design enables dynamic bit-width allocation to accommodate custom floating-point formats. The approach substantially reduces memory footprint and improves memory access efficiency. Experimental results demonstrate a 1.63× speedup over cuSPARSE in FP16 mode, while a custom format achieves FP32-level accuracy with performance surpassing FP16 cuSPARSE. Furthermore, a mixed-precision preconditioned conjugate gradient (PCG) solver attains a 2.09× acceleration.

This work addresses the challenge of automatically generating high-performance GPU tiled kernels from high-level tensor algebra expressions, thereby alleviating the burden of manual optimization. The authors propose an end-to-end compilation framework that integrates layer-wise lowering, expression rewriting, automated schedule search, reduction fusion, and tiling optimizations. For the first time, this framework automatically discovers high-efficiency kernels—comparable to FlashAttention-3—from the mathematical specification of attention operators, while introducing a novel scheduler that preserves program structural regularity. Evaluated on GH200 and RTX 5090 GPUs, the generated kernels achieve up to 23% and 42% higher throughput, respectively, and match or surpass hand-optimized cuDNN kernels across multiple long-sequence configurations.

This work addresses the underutilization of CPU resources in modern heterogeneous high-performance computing (HPC) systems when solving large-scale symmetric positive-definite linear systems using GPU-only approaches. Leveraging the SYCL programming model, the authors present the first heterogeneous implementations of the conjugate gradient (CG) method and Cholesky decomposition that operate across multi-vendor CPU-GPU platforms, including NVIDIA, AMD, and Intel architectures. Experimental results demonstrate that the heterogeneous CG solver achieves up to 32% speedup over GPU-only execution on large matrices, while the heterogeneous Cholesky decomposition attains a 29% acceleration. Furthermore, across diverse hardware vendors, the Cholesky solver consistently delivers at least a 12% performance improvement, significantly enhancing both computational efficiency and portability.