From 442b0878cebd7c2a0011f048de001c0b0c8c8fe2 Mon Sep 17 00:00:00 2001 From: Michel Schanen Date: Sat, 1 Aug 2026 21:30:00 -0500 Subject: [PATCH 1/2] Fix transposed result from sparse CSC matrix multiplication MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `dA * x` and `dA * B` returned Aᵀx and AᵀB for a oneSparseMatrixCSC. A CSC handle stores CSR(Aᵀ), so sparse_gemv! and sparse_gemm! apply flip_trans when calling into oneMKL and already present op(A) semantics to their callers. generic_matvecmul! and generic_matmatmul! applied flip_trans a second time, so 'N' reached oneMKL unchanged and it computed S*x = Aᵀx. Pass the transpose flag through instead. The existing sparse tests call sparse_gemv!/sparse_gemm! directly and never exercised `*` or `mul!` on a sparse matrix, which is why this went unnoticed; add a testset covering that path for both CSR and CSC. --- lib/mkl/interfaces.jl | 6 ++++-- test/onemkl.jl | 29 +++++++++++++++++++++++++++++ 2 files changed, 33 insertions(+), 2 deletions(-) diff --git a/lib/mkl/interfaces.jl b/lib/mkl/interfaces.jl index a7f43da4..77ea8728 100644 --- a/lib/mkl/interfaces.jl +++ b/lib/mkl/interfaces.jl @@ -18,7 +18,8 @@ function LinearAlgebra.generic_matvecmul!(C::oneVector{T}, tA::AbstractChar, A:: end function LinearAlgebra.generic_matvecmul!(C::oneVector{T}, tA::AbstractChar, A::oneSparseMatrixCSC{T}, B::oneVector{T}, alpha::Number, beta::Number) where {T <: BlasReal} - tA = tA in ('S', 's', 'H', 'h') ? 'T' : flip_trans(tA) + # sparse_gemv! already maps op(A) onto the transposed CSR handle, so tA is passed through + tA = tA in ('S', 's', 'H', 'h') ? 'N' : tA return sparse_gemv!(tA, alpha, A, B, beta, C) end @@ -29,7 +30,8 @@ function LinearAlgebra.generic_matmatmul!(C::oneMatrix{T}, tA, tB, A::oneSparseM end function LinearAlgebra.generic_matmatmul!(C::oneMatrix{T}, tA, tB, A::oneSparseMatrixCSC{T}, B::oneMatrix{T}, alpha::Number, beta::Number) where {T <: BlasReal} - tA = tA in ('S', 's', 'H', 'h') ? 'T' : flip_trans(tA) + # sparse_gemm! already maps op(A) onto the transposed CSR handle, so tA is passed through + tA = tA in ('S', 's', 'H', 'h') ? 'N' : tA tB = tB in ('S', 's', 'H', 'h') ? 'N' : tB return sparse_gemm!(tA, tB, alpha, A, B, beta, C) end diff --git a/test/onemkl.jl b/test/onemkl.jl index 31cd7b5f..4ec31840 100644 --- a/test/onemkl.jl +++ b/test/onemkl.jl @@ -1175,6 +1175,35 @@ end end end + # `*` and `mul!` reach oneMKL through generic_matvecmul!/generic_matmatmul!, which is + # a separate mapping from the sparse_gemv!/sparse_gemm! calls exercised above + @testset "sparse LinearAlgebra mul" begin + @testset "$SparseMatrix" for SparseMatrix in csr_csc_matrices + # the CSC methods are restricted to real element types + (SparseMatrix == oneSparseMatrixCSC && T <: Complex) && continue + + A = sprand(T, 10, 10, 0.5) + x = rand(T, 10) + y = rand(T, 10) + B = rand(T, 10, 2) + C = rand(T, 10, 2) + + dA = SparseMatrix(A) + dx = oneVector{T}(x) + dB = oneMatrix{T}(B) + + @test A * x ≈ collect(dA * dx) + @test A * B ≈ collect(dA * dB) + @test transpose(A) * x ≈ collect(transpose(dA) * dx) + @test transpose(A) * B ≈ collect(transpose(dA) * dB) + + alpha = rand(T) + beta = rand(T) + @test alpha * A * x + beta * y ≈ collect(mul!(oneVector{T}(y), dA, dx, alpha, beta)) + @test alpha * A * B + beta * C ≈ collect(mul!(oneMatrix{T}(C), dA, dB, alpha, beta)) + end + end + @testset "sparse symv" begin @testset "$SparseMatrix" for SparseMatrix in csr_csc_matrices @testset "uplo = $uplo" for uplo in ('L', 'U') From 1f20d191afe86378c06bafd97814620ac3613133 Mon Sep 17 00:00:00 2001 From: Michel Schanen Date: Sat, 1 Aug 2026 21:41:25 -0500 Subject: [PATCH 2/2] Support complex element types in sparse CSC multiplication The CSC generic_matvecmul!/generic_matmatmul! methods were constrained to BlasReal, so `dA * x` on a complex oneSparseMatrixCSC never reached them and fell through to the AbstractArray fallback, failing with CanonicalIndexError: getindex not defined for oneSparseMatrixCSC{ComplexF32, Int64} The underlying sparse_gemv!/sparse_gemm! wrappers already carry complex CSC specializations that handle the conjugation identities for the 'C' case, so widening the constraint to BlasFloat, matching the CSR methods, is all that is needed. Extend the multiplication testset over identity/transpose/adjoint for every supported element type, which covers the complex-only distinction between transpose and adjoint. --- lib/mkl/interfaces.jl | 4 ++-- test/onemkl.jl | 11 ++++------- 2 files changed, 6 insertions(+), 9 deletions(-) diff --git a/lib/mkl/interfaces.jl b/lib/mkl/interfaces.jl index 77ea8728..a18f74ea 100644 --- a/lib/mkl/interfaces.jl +++ b/lib/mkl/interfaces.jl @@ -17,7 +17,7 @@ function LinearAlgebra.generic_matvecmul!(C::oneVector{T}, tA::AbstractChar, A:: return sparse_gemv!(tA, alpha, A, B, beta, C) end -function LinearAlgebra.generic_matvecmul!(C::oneVector{T}, tA::AbstractChar, A::oneSparseMatrixCSC{T}, B::oneVector{T}, alpha::Number, beta::Number) where {T <: BlasReal} +function LinearAlgebra.generic_matvecmul!(C::oneVector{T}, tA::AbstractChar, A::oneSparseMatrixCSC{T}, B::oneVector{T}, alpha::Number, beta::Number) where {T <: BlasFloat} # sparse_gemv! already maps op(A) onto the transposed CSR handle, so tA is passed through tA = tA in ('S', 's', 'H', 'h') ? 'N' : tA return sparse_gemv!(tA, alpha, A, B, beta, C) @@ -29,7 +29,7 @@ function LinearAlgebra.generic_matmatmul!(C::oneMatrix{T}, tA, tB, A::oneSparseM return sparse_gemm!(tA, tB, alpha, A, B, beta, C) end -function LinearAlgebra.generic_matmatmul!(C::oneMatrix{T}, tA, tB, A::oneSparseMatrixCSC{T}, B::oneMatrix{T}, alpha::Number, beta::Number) where {T <: BlasReal} +function LinearAlgebra.generic_matmatmul!(C::oneMatrix{T}, tA, tB, A::oneSparseMatrixCSC{T}, B::oneMatrix{T}, alpha::Number, beta::Number) where {T <: BlasFloat} # sparse_gemm! already maps op(A) onto the transposed CSR handle, so tA is passed through tA = tA in ('S', 's', 'H', 'h') ? 'N' : tA tB = tB in ('S', 's', 'H', 'h') ? 'N' : tB diff --git a/test/onemkl.jl b/test/onemkl.jl index 4ec31840..f0383962 100644 --- a/test/onemkl.jl +++ b/test/onemkl.jl @@ -1179,9 +1179,6 @@ end # a separate mapping from the sparse_gemv!/sparse_gemm! calls exercised above @testset "sparse LinearAlgebra mul" begin @testset "$SparseMatrix" for SparseMatrix in csr_csc_matrices - # the CSC methods are restricted to real element types - (SparseMatrix == oneSparseMatrixCSC && T <: Complex) && continue - A = sprand(T, 10, 10, 0.5) x = rand(T, 10) y = rand(T, 10) @@ -1192,10 +1189,10 @@ end dx = oneVector{T}(x) dB = oneMatrix{T}(B) - @test A * x ≈ collect(dA * dx) - @test A * B ≈ collect(dA * dB) - @test transpose(A) * x ≈ collect(transpose(dA) * dx) - @test transpose(A) * B ≈ collect(transpose(dA) * dB) + @testset "opa = $(nameof(opa))" for opa in (identity, transpose, adjoint) + @test opa(A) * x ≈ collect(opa(dA) * dx) + @test opa(A) * B ≈ collect(opa(dA) * dB) + end alpha = rand(T) beta = rand(T)