From e7db0150e44acafed78059600a73691005fa3aca Mon Sep 17 00:00:00 2001 From: Anton Pozharskiy Date: Fri, 21 Aug 2026 21:35:27 +0200 Subject: [PATCH 1/9] fix bug in transpose, add symmetric multiply --- src/level2.jl | 65 +++++++++++++++++++++++++++++++++++++++++++--- src/vecs.jl | 1 + test/mat/level2.jl | 30 +++++++++++++++++++++ 3 files changed, 92 insertions(+), 4 deletions(-) diff --git a/src/level2.jl b/src/level2.jl index df279e4..7486604 100644 --- a/src/level2.jl +++ b/src/level2.jl @@ -11,12 +11,14 @@ for (Mat, Vec, flag) in [ blasfeo_gemm_nt = Symbol(:blasfeo_, flag, :gemm_nt) blasfeo_gemm_tn = Symbol(:blasfeo_, flag, :gemm_tn) blasfeo_gemm_tt = Symbol(:blasfeo_, flag, :gemm_tt) + blasfeo_symv_l = Symbol(:blasfeo_, flag, :symv_l) + blasfeo_symv_u = Symbol(:blasfeo_, flag, :symv_u) @eval function Base.:*(A::$Mat, x::$Vec) @boundscheck begin size(A)[2] == length(x) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) end - z = similar(x) + z = similar(x, size(A,1)) $blasfeo_gemv_n( size(A,1), size(A,2), 1.0, A, 0, 0, @@ -29,12 +31,12 @@ for (Mat, Vec, flag) in [ @eval function Base.:*(A::Transpose{eltype($Mat), $Mat}, x::$Vec) @boundscheck begin - size(A.parent)[1] == length(x) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) + size(A,2) == length(x) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) end - z = similar(x) + z = similar(x,size(A,1)) $blasfeo_gemv_t( - size(A,1), size(A,2), + size(A,2), size(A,1), 1.0, A, 0, 0, x, 0, 0.0, x, 0, @@ -73,4 +75,59 @@ for (Mat, Vec, flag) in [ ) return Y end + + + @eval function Base.:*(A::Symmetric{eltype($Mat), $Mat}, x::$Vec) + @boundscheck begin + size(A,1) == length(x) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) + end + + z = similar(x) + if A.uplo == 'L' + $blasfeo_symv_l( + size(A,1), + 1.0, A.data, 0, 0, + x, 0, + 0.0, x, 0, + z, 0, + ) + else + $blasfeo_symv_u( + size(A,1), + 1.0, A.data, 0, 0, + x, 0, + 0.0, x, 0, + z, 0, + ) + end + return z # A^T*x + end + + @eval function LinearAlgebra.mul!(Y::$Vec, A::Symmetric{eltype($Mat), $Mat}, B::$Vec) + @boundscheck begin + size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) + size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) + end + + if A.uplo == 'L' + $blasfeo_symv_l( + size(A,1), + 1.0, A.data, 0, 0, + B, 0, + 0.0, B, 0, + Y, 0, + ) + else + $blasfeo_symv_u( + size(A,1), + 1.0, A.data, 0, 0, + B, 0, + 0.0, B, 0, + Y, 0, + ) + end + return Y + end + + end diff --git a/src/vecs.jl b/src/vecs.jl index 1704d00..584d37f 100644 --- a/src/vecs.jl +++ b/src/vecs.jl @@ -81,6 +81,7 @@ for (type,flag) in [ # similar @eval Base.similar(A::$type) = $type(length(A)) @eval Base.similar(A::$type, dims::Dims{1}) = $type(dims...) + @eval Base.similar(A::$type, dims::Integer) = $type(dims) @eval Base.similar(A::$type, ::eltype($type)) = $type(length(A)) @eval Base.similar(A::$type, ::eltype($type), dims::Dims{1}) = $type(dims...) diff --git a/test/mat/level2.jl b/test/mat/level2.jl index d1e275b..e14471e 100644 --- a/test/mat/level2.jl +++ b/test/mat/level2.jl @@ -26,5 +26,35 @@ # test mul! @test A*a ≈ mul!(y_blasfeo, A_blasfeo, a_blasfeo) @test transpose(A)*a ≈ mul!(y_blasfeo, transpose(A_blasfeo), a_blasfeo) + + # test asymmetric + A = rand(eltype(VEC), 100, 50) + a = rand(eltype(VEC), 50) + y = rand(eltype(VEC), 100) + A_blasfeo = MAT(A) + a_blasfeo = VEC(a) + y_blasfeo = VEC(y) + + @test A*a ≈ A_blasfeo*a_blasfeo + @test isa(A_blasfeo*a_blasfeo, VEC) + @test transpose(A)*y ≈ transpose(A_blasfeo)*y_blasfeo + @test isa(transpose(A_blasfeo)*y_blasfeo, VEC) + + # test symmetric + A = rand(eltype(VEC), 100, 100) + a = rand(eltype(VEC), 100) + y = rand(eltype(VEC), 100) + A_blasfeo = MAT(A) + a_blasfeo = VEC(a) + y_blasfeo = VEC(y) + A_L = Symmetric(A,:L) + A_L_blasfeo = Symmetric(A_blasfeo,:L) + A_U = Symmetric(A,:U) + A_U_blasfeo = Symmetric(A_blasfeo,:U) + + @test A_L*a ≈ A_L_blasfeo*a_blasfeo + @test isa(A_L_blasfeo*a_blasfeo, VEC) + @test A_U*a ≈ A_U_blasfeo*a_blasfeo + @test isa(A_U_blasfeo*a_blasfeo, VEC) end end From 5f007562853ee5d58756f988389a86a0fdbfbe78 Mon Sep 17 00:00:00 2001 From: Anton Pozharskiy Date: Fri, 21 Aug 2026 22:51:35 +0200 Subject: [PATCH 2/9] triangular matrices --- src/level2.jl | 135 +++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 133 insertions(+), 2 deletions(-) diff --git a/src/level2.jl b/src/level2.jl index 7486604..ab88d55 100644 --- a/src/level2.jl +++ b/src/level2.jl @@ -13,6 +13,10 @@ for (Mat, Vec, flag) in [ blasfeo_gemm_tt = Symbol(:blasfeo_, flag, :gemm_tt) blasfeo_symv_l = Symbol(:blasfeo_, flag, :symv_l) blasfeo_symv_u = Symbol(:blasfeo_, flag, :symv_u) + blasfeo_trmv_lnn = Symbol(:blasfeo_, flag, :trmv_lnn) + blasfeo_trmv_lnu = Symbol(:blasfeo_, flag, :trmv_lnu) + blasfeo_trmv_unn = Symbol(:blasfeo_, flag, :trmv_unn) + blasfeo_trmv_unu = Symbol(:blasfeo_, flag, :trmv_unu) @eval function Base.:*(A::$Mat, x::$Vec) @boundscheck begin @@ -100,7 +104,7 @@ for (Mat, Vec, flag) in [ z, 0, ) end - return z # A^T*x + return z end @eval function LinearAlgebra.mul!(Y::$Vec, A::Symmetric{eltype($Mat), $Mat}, B::$Vec) @@ -129,5 +133,132 @@ for (Mat, Vec, flag) in [ return Y end - + + @eval function Base.:*(A::LowerTriangular{eltype($Mat), $Mat}, x::$Vec) + @boundscheck begin + size(A,1) == length(x) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) + end + + z = similar(x) + $blasfeo_trmv_lnn( + size(A,1), + 1.0, A.data, 0, 0, + x, 0, + 0.0, x, 0, + z, 0, + ) + return z # A^T*x + end + + @eval function LinearAlgebra.mul!(Y::$Vec, A::LowerTriangular{eltype($Mat), $Mat}, B::$Vec) + @boundscheck begin + size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) + size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) + end + + $blasfeo_trmv_lnn( + size(A,1), + 1.0, A.data, 0, 0, + B, 0, + 0.0, B, 0, + Y, 0, + ) + return Y + end + + @eval function Base.:*(A::UpperTriangular{eltype($Mat), $Mat}, x::$Vec) + @boundscheck begin + size(A,1) == length(x) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) + end + + z = similar(x) + $blasfeo_trmv_unn( + size(A,1), + 1.0, A.data, 0, 0, + x, 0, + 0.0, x, 0, + z, 0, + ) + return z # A^T*x + end + + @eval function LinearAlgebra.mul!(Y::$Vec, A::UpperTriangular{eltype($Mat), $Mat}, B::$Vec) + @boundscheck begin + size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) + size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) + end + + $blasfeo_trmv_unn( + size(A,1), + 1.0, A.data, 0, 0, + B, 0, + 0.0, B, 0, + Y, 0, + ) + return Y + end + + @eval function Base.:*(A::UnitLowerTriangular{eltype($Mat), $Mat}, x::$Vec) + @boundscheck begin + size(A,1) == length(x) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) + end + + z = similar(x) + $blasfeo_trmv_lnu( + size(A,1), + 1.0, A.data, 0, 0, + x, 0, + 0.0, x, 0, + z, 0, + ) + return z # A^T*x + end + + @eval function LinearAlgebra.mul!(Y::$Vec, A::UnitLowerTriangular{eltype($Mat), $Mat}, B::$Vec) + @boundscheck begin + size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) + size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) + end + + $blasfeo_trmv_lnu( + size(A,1), + 1.0, A.data, 0, 0, + B, 0, + 0.0, B, 0, + Y, 0, + ) + return Y + end + + @eval function Base.:*(A::UnitUpperTriangular{eltype($Mat), $Mat}, x::$Vec) + @boundscheck begin + size(A,1) == length(x) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) + end + + z = similar(x) + $blasfeo_trmv_unu( + size(A,1), + 1.0, A.data, 0, 0, + x, 0, + 0.0, x, 0, + z, 0, + ) + return z # A^T*x + end + + @eval function LinearAlgebra.mul!(Y::$Vec, A::UnitUpperTriangular{eltype($Mat), $Mat}, B::$Vec) + @boundscheck begin + size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) + size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) + end + + $blasfeo_trmv_unu( + size(A,1), + 1.0, A.data, 0, 0, + B, 0, + 0.0, B, 0, + Y, 0, + ) + return Y + end end From 9360fc972cab3059faec005323c3f97e96b538b1 Mon Sep 17 00:00:00 2001 From: Anton Pozharskiy Date: Sat, 22 Aug 2026 12:19:35 +0200 Subject: [PATCH 3/9] non-transposed trmv tests+workarounds for unimplemented versions --- src/level2.jl | 222 ++++++++++++++------------------------------- test/mat/level2.jl | 90 +++++++++++------- 2 files changed, 124 insertions(+), 188 deletions(-) diff --git a/src/level2.jl b/src/level2.jl index ab88d55..02be8ee 100644 --- a/src/level2.jl +++ b/src/level2.jl @@ -1,7 +1,7 @@ -for (Mat, Vec, flag) in [ - (:BlasfeoDmat, :BlasfeoDvec, :d), - (:BlasfeoSmat, :BlasfeoSvec, :s), +for (El, Mat, Vec, flag) in [ + (:Cdouble, :BlasfeoDmat, :BlasfeoDvec, :d), + (:Cfloat, :BlasfeoSmat, :BlasfeoSvec, :s), ] # Overload matrix-vector multiplication @@ -19,34 +19,13 @@ for (Mat, Vec, flag) in [ blasfeo_trmv_unu = Symbol(:blasfeo_, flag, :trmv_unu) @eval function Base.:*(A::$Mat, x::$Vec) - @boundscheck begin - size(A)[2] == length(x) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) - end z = similar(x, size(A,1)) - $blasfeo_gemv_n( - size(A,1), size(A,2), - 1.0, A, 0, 0, - x, 0, - 0.0, x, 0, - z, 0, - ) - return z # A*x + return mul!(z,A,x) end - @eval function Base.:*(A::Transpose{eltype($Mat), $Mat}, x::$Vec) - @boundscheck begin - size(A,2) == length(x) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) - end - - z = similar(x,size(A,1)) - $blasfeo_gemv_t( - size(A,2), size(A,1), - 1.0, A, 0, 0, - x, 0, - 0.0, x, 0, - z, 0, - ) - return z # A^T*x + @eval function Base.:*(A::Transpose{$El, $Mat}, x::$Vec) + z = similar(x, size(A,1)) + return mul!(z,A,x) # A^T*x end @eval function LinearAlgebra.mul!(Y::$Vec, A::$Mat, B::$Vec) @@ -64,14 +43,14 @@ for (Mat, Vec, flag) in [ return Y end - @eval function LinearAlgebra.mul!(Y::$Vec, A::Transpose{eltype($Mat), $Mat}, B::$Vec) + @eval function LinearAlgebra.mul!(Y::$Vec, A::Transpose{$El, $Mat}, B::$Vec) @boundscheck begin size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) end $blasfeo_gemv_t( - size(A,1), size(A,2), + size(A,2), size(A,1), 1.0, A, 0, 0, B, 0, 0.0, B, 0, @@ -80,34 +59,12 @@ for (Mat, Vec, flag) in [ return Y end - - @eval function Base.:*(A::Symmetric{eltype($Mat), $Mat}, x::$Vec) - @boundscheck begin - size(A,1) == length(x) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) - end - + @eval function Base.:*(A::Symmetric{$El, $Mat}, x::$Vec) z = similar(x) - if A.uplo == 'L' - $blasfeo_symv_l( - size(A,1), - 1.0, A.data, 0, 0, - x, 0, - 0.0, x, 0, - z, 0, - ) - else - $blasfeo_symv_u( - size(A,1), - 1.0, A.data, 0, 0, - x, 0, - 0.0, x, 0, - z, 0, - ) - end - return z + return mul!(z,A,x) end - @eval function LinearAlgebra.mul!(Y::$Vec, A::Symmetric{eltype($Mat), $Mat}, B::$Vec) + @eval function LinearAlgebra.mul!(Y::$Vec, A::Symmetric{$El, $Mat}, B::$Vec) @boundscheck begin size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) @@ -134,23 +91,12 @@ for (Mat, Vec, flag) in [ end - @eval function Base.:*(A::LowerTriangular{eltype($Mat), $Mat}, x::$Vec) - @boundscheck begin - size(A,1) == length(x) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) - end - + @eval function Base.:*(A::LowerTriangular{$El, $Mat}, x::$Vec) z = similar(x) - $blasfeo_trmv_lnn( - size(A,1), - 1.0, A.data, 0, 0, - x, 0, - 0.0, x, 0, - z, 0, - ) - return z # A^T*x + return mul!(z,A,x) end - @eval function LinearAlgebra.mul!(Y::$Vec, A::LowerTriangular{eltype($Mat), $Mat}, B::$Vec) + @eval function LinearAlgebra.mul!(Y::$Vec, A::LowerTriangular{$El, $Mat}, B::$Vec) @boundscheck begin size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) @@ -158,107 +104,73 @@ for (Mat, Vec, flag) in [ $blasfeo_trmv_lnn( size(A,1), - 1.0, A.data, 0, 0, + A.data, 0, 0, B, 0, - 0.0, B, 0, Y, 0, ) return Y end - @eval function Base.:*(A::UpperTriangular{eltype($Mat), $Mat}, x::$Vec) - @boundscheck begin - size(A,1) == length(x) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) - end - - z = similar(x) - $blasfeo_trmv_unn( - size(A,1), - 1.0, A.data, 0, 0, - x, 0, - 0.0, x, 0, - z, 0, - ) - return z # A^T*x - end - - @eval function LinearAlgebra.mul!(Y::$Vec, A::UpperTriangular{eltype($Mat), $Mat}, B::$Vec) - @boundscheck begin - size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) - size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) + if El == :Cdouble # TODO(@anton) blasfeo_strmv_unn and blasfeo_strmv_lnu are unimplemented :( + @eval function Base.:*(A::UpperTriangular{$El, $Mat}, x::$Vec) + z = similar(x) + return mul!(z,A,x) end - $blasfeo_trmv_unn( - size(A,1), - 1.0, A.data, 0, 0, - B, 0, - 0.0, B, 0, - Y, 0, - ) - return Y - end + @eval function LinearAlgebra.mul!(Y::$Vec, A::UpperTriangular{$El, $Mat}, B::$Vec) + @boundscheck begin + size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) + size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) + end - @eval function Base.:*(A::UnitLowerTriangular{eltype($Mat), $Mat}, x::$Vec) - @boundscheck begin - size(A,1) == length(x) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) + $blasfeo_trmv_unn( + size(A,1), + A.data, 0, 0, + B, 0, + Y, 0, + ) + return Y end - z = similar(x) - $blasfeo_trmv_lnu( - size(A,1), - 1.0, A.data, 0, 0, - x, 0, - 0.0, x, 0, - z, 0, - ) - return z # A^T*x - end - - @eval function LinearAlgebra.mul!(Y::$Vec, A::UnitLowerTriangular{eltype($Mat), $Mat}, B::$Vec) - @boundscheck begin - size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) - size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) + @eval function Base.:*(A::UnitLowerTriangular{$El, $Mat}, x::$Vec) + z = similar(x) + return mul!(z,A,x) end - $blasfeo_trmv_lnu( - size(A,1), - 1.0, A.data, 0, 0, - B, 0, - 0.0, B, 0, - Y, 0, - ) - return Y - end + @eval function LinearAlgebra.mul!(Y::$Vec, A::UnitLowerTriangular{$El, $Mat}, B::$Vec) + @boundscheck begin + size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) + size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) + end - @eval function Base.:*(A::UnitUpperTriangular{eltype($Mat), $Mat}, x::$Vec) - @boundscheck begin - size(A,1) == length(x) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) + $blasfeo_trmv_lnu( + size(A,1),A.data, 0, 0, + B, 0, + Y, 0, + ) + return Y end - - z = similar(x) - $blasfeo_trmv_unu( - size(A,1), - 1.0, A.data, 0, 0, - x, 0, - 0.0, x, 0, - z, 0, - ) - return z # A^T*x end - @eval function LinearAlgebra.mul!(Y::$Vec, A::UnitUpperTriangular{eltype($Mat), $Mat}, B::$Vec) - @boundscheck begin - size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) - size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) - end - - $blasfeo_trmv_unu( - size(A,1), - 1.0, A.data, 0, 0, - B, 0, - 0.0, B, 0, - Y, 0, - ) - return Y - end + # need to implement trmv_unu + # @eval function Base.:*(A::UnitUpperTriangular{$El, $Mat}, x::$Vec) + # z = similar(x) + # return mul!(z,A,x) + # end + + # @eval function LinearAlgebra.mul!(Y::$Vec, A::UnitUpperTriangular{$El, $Mat}, B::$Vec) + # @boundscheck begin + # size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) + # size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) + # end + + # $blasfeo_trmv_unu( + # size(A,1), + # 1.0, A.data, 0, 0, + # B, 0, + # 0.0, B, 0, + # Y, 0, + # ) + # return Y + # end end diff --git a/test/mat/level2.jl b/test/mat/level2.jl index e14471e..cf62fe4 100644 --- a/test/mat/level2.jl +++ b/test/mat/level2.jl @@ -1,15 +1,26 @@ @testset "Level 2 BLAS Operations" begin @testset for (MAT,VEC) in ((BlasfeoDmat,BlasfeoDvec), (BlasfeoSmat,BlasfeoSvec)) - A = rand(eltype(VEC), 100, 100) - B = rand(eltype(VEC), 95, 95) - a = rand(eltype(VEC), 100) - y = rand(eltype(VEC), 100) - b = rand(eltype(VEC), 90) + A_orig = rand(eltype(VEC), 100, 100) + B_orig = rand(eltype(VEC), 95, 95) + C_orig = rand(eltype(VEC), 100, 50) + a_orig = rand(eltype(VEC), 100) + y_orig = rand(eltype(VEC), 100) + b_orig = rand(eltype(VEC), 90) + c_orig = rand(eltype(VEC), 50) + A = copy(A_orig) + B = copy(B_orig) + C = copy(C_orig) + a = copy(a_orig) + y = copy(y_orig) + b = copy(b_orig) + c = copy(c_orig) A_blasfeo = MAT(A) B_blasfeo = MAT(B) + C_blasfeo = MAT(C) a_blasfeo = VEC(a) y_blasfeo = VEC(y) b_blasfeo = VEC(b) + c_blasfeo = VEC(c) @test_throws DimensionMismatch A_blasfeo*b_blasfeo @test_throws DimensionMismatch B_blasfeo*a_blasfeo @@ -24,37 +35,50 @@ @test isa(transpose(A_blasfeo)*a_blasfeo, VEC) # test mul! - @test A*a ≈ mul!(y_blasfeo, A_blasfeo, a_blasfeo) - @test transpose(A)*a ≈ mul!(y_blasfeo, transpose(A_blasfeo), a_blasfeo) + @test mul!(y,A,a) ≈ mul!(y_blasfeo, A_blasfeo, a_blasfeo) + @test mul!(y,transpose(A),a) ≈ mul!(y_blasfeo, transpose(A_blasfeo), a_blasfeo) # test asymmetric - A = rand(eltype(VEC), 100, 50) - a = rand(eltype(VEC), 50) - y = rand(eltype(VEC), 100) - A_blasfeo = MAT(A) - a_blasfeo = VEC(a) - y_blasfeo = VEC(y) - - @test A*a ≈ A_blasfeo*a_blasfeo - @test isa(A_blasfeo*a_blasfeo, VEC) - @test transpose(A)*y ≈ transpose(A_blasfeo)*y_blasfeo - @test isa(transpose(A_blasfeo)*y_blasfeo, VEC) + @test C*c ≈ C_blasfeo*c_blasfeo + @test isa(C_blasfeo*c_blasfeo, VEC) + @test transpose(C)*y ≈ transpose(C_blasfeo)*y_blasfeo + @test isa(transpose(C_blasfeo)*y_blasfeo, VEC) # test symmetric - A = rand(eltype(VEC), 100, 100) - a = rand(eltype(VEC), 100) - y = rand(eltype(VEC), 100) - A_blasfeo = MAT(A) - a_blasfeo = VEC(a) - y_blasfeo = VEC(y) - A_L = Symmetric(A,:L) - A_L_blasfeo = Symmetric(A_blasfeo,:L) - A_U = Symmetric(A,:U) - A_U_blasfeo = Symmetric(A_blasfeo,:U) - - @test A_L*a ≈ A_L_blasfeo*a_blasfeo - @test isa(A_L_blasfeo*a_blasfeo, VEC) - @test A_U*a ≈ A_U_blasfeo*a_blasfeo - @test isa(A_U_blasfeo*a_blasfeo, VEC) + A_sym_L = Symmetric(A,:L) + A_sym_L_blasfeo = Symmetric(A_blasfeo,:L) + A_sym_U = Symmetric(A,:U) + A_sym_U_blasfeo = Symmetric(A_blasfeo,:U) + + @test A_sym_L*a ≈ A_sym_L_blasfeo*a_blasfeo + @test isa(A_sym_L_blasfeo*a_blasfeo, VEC) + @test A_sym_U*a ≈ A_sym_U_blasfeo*a_blasfeo + @test isa(A_sym_U_blasfeo*a_blasfeo, VEC) + + # test triangular + A_lnn = LowerTriangular(A) + A_lnn_blasfeo = LowerTriangular(A_blasfeo) + A_unn = UpperTriangular(A) + A_unn_blasfeo = UpperTriangular(A_blasfeo) + A_lnu = UnitLowerTriangular(A) + A_lnu_blasfeo = UnitLowerTriangular(A_blasfeo) + A_unu = UnitUpperTriangular(A) + A_unu_blasfeo = UnitUpperTriangular(A_blasfeo) + + # test lnn + @test A_lnn*a ≈ A_lnn_blasfeo*a_blasfeo + @test isa(A_lnn_blasfeo*a_blasfeo, VEC) + if MAT == BlasfeoDmat # TODO(@anton) not implemented upstream + # test unn + @test A_unn*a ≈ A_unn_blasfeo*a_blasfeo + @test isa(A_unn_blasfeo*a_blasfeo, VEC) + # test lnu + @test A_lnu*a ≈ A_lnu_blasfeo*a_blasfeo + @test isa(A_lnu_blasfeo*a_blasfeo, VEC) + # test unu + # TODO(@anton) not implemented upstream + #@test A_unu*a ≈ A_unu_blasfeo*a_blasfeo + #@test isa(A_unu_blasfeo*a_blasfeo, VEC) + end end end From 1eab1a24152b9f37927242d0f37bd7ac309ca6b0 Mon Sep 17 00:00:00 2001 From: Anton Pozharskiy Date: Sun, 23 Aug 2026 15:40:11 +0200 Subject: [PATCH 4/9] separating level 2 implementations and working part of triangular multiplies working --- src/level2.jl | 179 +-------------------------------------------- src/level2/gemv.jl | 48 ++++++++++++ src/level2/symv.jl | 39 ++++++++++ src/level2/trmv.jl | 155 +++++++++++++++++++++++++++++++++++++++ test/mat/level2.jl | 164 +++++++++++++++++++++++------------------ 5 files changed, 338 insertions(+), 247 deletions(-) create mode 100644 src/level2/gemv.jl create mode 100644 src/level2/symv.jl create mode 100644 src/level2/trmv.jl diff --git a/src/level2.jl b/src/level2.jl index 02be8ee..b824d08 100644 --- a/src/level2.jl +++ b/src/level2.jl @@ -1,176 +1,3 @@ - -for (El, Mat, Vec, flag) in [ - (:Cdouble, :BlasfeoDmat, :BlasfeoDvec, :d), - (:Cfloat, :BlasfeoSmat, :BlasfeoSvec, :s), - ] - # Overload matrix-vector multiplication - - blasfeo_gemv_n = Symbol(:blasfeo_, flag, :gemv_n) - blasfeo_gemv_t = Symbol(:blasfeo_, flag, :gemv_t) - blasfeo_gemm_nn = Symbol(:blasfeo_, flag, :gemm_nn) - blasfeo_gemm_nt = Symbol(:blasfeo_, flag, :gemm_nt) - blasfeo_gemm_tn = Symbol(:blasfeo_, flag, :gemm_tn) - blasfeo_gemm_tt = Symbol(:blasfeo_, flag, :gemm_tt) - blasfeo_symv_l = Symbol(:blasfeo_, flag, :symv_l) - blasfeo_symv_u = Symbol(:blasfeo_, flag, :symv_u) - blasfeo_trmv_lnn = Symbol(:blasfeo_, flag, :trmv_lnn) - blasfeo_trmv_lnu = Symbol(:blasfeo_, flag, :trmv_lnu) - blasfeo_trmv_unn = Symbol(:blasfeo_, flag, :trmv_unn) - blasfeo_trmv_unu = Symbol(:blasfeo_, flag, :trmv_unu) - - @eval function Base.:*(A::$Mat, x::$Vec) - z = similar(x, size(A,1)) - return mul!(z,A,x) - end - - @eval function Base.:*(A::Transpose{$El, $Mat}, x::$Vec) - z = similar(x, size(A,1)) - return mul!(z,A,x) # A^T*x - end - - @eval function LinearAlgebra.mul!(Y::$Vec, A::$Mat, B::$Vec) - @boundscheck begin - size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) - size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) - end - $blasfeo_gemv_n( - size(A,1), size(A,2), - 1.0, A, 0, 0, - B, 0, - 0.0, B, 0, - Y, 0, - ) - return Y - end - - @eval function LinearAlgebra.mul!(Y::$Vec, A::Transpose{$El, $Mat}, B::$Vec) - @boundscheck begin - size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) - size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) - end - - $blasfeo_gemv_t( - size(A,2), size(A,1), - 1.0, A, 0, 0, - B, 0, - 0.0, B, 0, - Y, 0, - ) - return Y - end - - @eval function Base.:*(A::Symmetric{$El, $Mat}, x::$Vec) - z = similar(x) - return mul!(z,A,x) - end - - @eval function LinearAlgebra.mul!(Y::$Vec, A::Symmetric{$El, $Mat}, B::$Vec) - @boundscheck begin - size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) - size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) - end - - if A.uplo == 'L' - $blasfeo_symv_l( - size(A,1), - 1.0, A.data, 0, 0, - B, 0, - 0.0, B, 0, - Y, 0, - ) - else - $blasfeo_symv_u( - size(A,1), - 1.0, A.data, 0, 0, - B, 0, - 0.0, B, 0, - Y, 0, - ) - end - return Y - end - - - @eval function Base.:*(A::LowerTriangular{$El, $Mat}, x::$Vec) - z = similar(x) - return mul!(z,A,x) - end - - @eval function LinearAlgebra.mul!(Y::$Vec, A::LowerTriangular{$El, $Mat}, B::$Vec) - @boundscheck begin - size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) - size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) - end - - $blasfeo_trmv_lnn( - size(A,1), - A.data, 0, 0, - B, 0, - Y, 0, - ) - return Y - end - - if El == :Cdouble # TODO(@anton) blasfeo_strmv_unn and blasfeo_strmv_lnu are unimplemented :( - @eval function Base.:*(A::UpperTriangular{$El, $Mat}, x::$Vec) - z = similar(x) - return mul!(z,A,x) - end - - @eval function LinearAlgebra.mul!(Y::$Vec, A::UpperTriangular{$El, $Mat}, B::$Vec) - @boundscheck begin - size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) - size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) - end - - $blasfeo_trmv_unn( - size(A,1), - A.data, 0, 0, - B, 0, - Y, 0, - ) - return Y - end - - @eval function Base.:*(A::UnitLowerTriangular{$El, $Mat}, x::$Vec) - z = similar(x) - return mul!(z,A,x) - end - - @eval function LinearAlgebra.mul!(Y::$Vec, A::UnitLowerTriangular{$El, $Mat}, B::$Vec) - @boundscheck begin - size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) - size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) - end - - $blasfeo_trmv_lnu( - size(A,1),A.data, 0, 0, - B, 0, - Y, 0, - ) - return Y - end - end - - # need to implement trmv_unu - # @eval function Base.:*(A::UnitUpperTriangular{$El, $Mat}, x::$Vec) - # z = similar(x) - # return mul!(z,A,x) - # end - - # @eval function LinearAlgebra.mul!(Y::$Vec, A::UnitUpperTriangular{$El, $Mat}, B::$Vec) - # @boundscheck begin - # size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) - # size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) - # end - - # $blasfeo_trmv_unu( - # size(A,1), - # 1.0, A.data, 0, 0, - # B, 0, - # 0.0, B, 0, - # Y, 0, - # ) - # return Y - # end -end +include("level2/gemv.jl") +include("level2/symv.jl") +include("level2/trmv.jl") diff --git a/src/level2/gemv.jl b/src/level2/gemv.jl new file mode 100644 index 0000000..2f8137c --- /dev/null +++ b/src/level2/gemv.jl @@ -0,0 +1,48 @@ +for (El, Mat, Vec, flag) in [ + (:Cdouble, :BlasfeoDmat, :BlasfeoDvec, :d), + (:Cfloat, :BlasfeoSmat, :BlasfeoSvec, :s), + ] + blasfeo_gemv_n = Symbol(:blasfeo_, flag, :gemv_n) + blasfeo_gemv_t = Symbol(:blasfeo_, flag, :gemv_t) + + @eval function Base.:*(A::$Mat, x::$Vec) + z = similar(x, size(A,1)) + return mul!(z,A,x) + end + + @eval function Base.:*(A::Transpose{$El, $Mat}, x::$Vec) + z = similar(x, size(A,1)) + return mul!(z,A,x) # A^T*x + end + + @eval function LinearAlgebra.mul!(Y::$Vec, A::$Mat, B::$Vec) + @boundscheck begin + size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) + size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) + end + $blasfeo_gemv_n( + size(A,1), size(A,2), + 1.0, A, 0, 0, + B, 0, + 0.0, B, 0, + Y, 0, + ) + return Y + end + + @eval function LinearAlgebra.mul!(Y::$Vec, A::Transpose{$El, $Mat}, B::$Vec) + @boundscheck begin + size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) + size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) + end + $blasfeo_gemv_t( + size(A,2), size(A,1), + 1.0, A, 0, 0, + B, 0, + 0.0, B, 0, + Y, 0, + ) + return Y + end + +end diff --git a/src/level2/symv.jl b/src/level2/symv.jl new file mode 100644 index 0000000..e9fb186 --- /dev/null +++ b/src/level2/symv.jl @@ -0,0 +1,39 @@ +for (El, Mat, Vec, flag) in [ + (:Cdouble, :BlasfeoDmat, :BlasfeoDvec, :d), + (:Cfloat, :BlasfeoSmat, :BlasfeoSvec, :s), + ] + # Overload matrix-vector multiplication + blasfeo_symv_l = Symbol(:blasfeo_, flag, :symv_l) + blasfeo_symv_u = Symbol(:blasfeo_, flag, :symv_u) + + @eval function Base.:*(A::Symmetric{$El, $Mat}, x::$Vec) + z = similar(x) + return mul!(z,A,x) + end + + @eval function LinearAlgebra.mul!(Y::$Vec, A::Symmetric{$El, $Mat}, B::$Vec) + @boundscheck begin + size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) + size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) + end + + if A.uplo == 'L' + $blasfeo_symv_l( + size(A,1), + 1.0, A.data, 0, 0, + B, 0, + 0.0, B, 0, + Y, 0, + ) + else + $blasfeo_symv_u( + size(A,1), + 1.0, A.data, 0, 0, + B, 0, + 0.0, B, 0, + Y, 0, + ) + end + return Y + end +end diff --git a/src/level2/trmv.jl b/src/level2/trmv.jl new file mode 100644 index 0000000..461dd89 --- /dev/null +++ b/src/level2/trmv.jl @@ -0,0 +1,155 @@ +for (El, Mat, Vec, flag) in [ + (:Cdouble, :BlasfeoDmat, :BlasfeoDvec, :d), + (:Cfloat, :BlasfeoSmat, :BlasfeoSvec, :s), + ] + + blasfeo_trmv_lnn = Symbol(:blasfeo_, flag, :trmv_lnn) + blasfeo_trmv_ltn = Symbol(:blasfeo_, flag, :trmv_ltn) + blasfeo_trmv_lnu = Symbol(:blasfeo_, flag, :trmv_lnu) + blasfeo_trmv_ltu = Symbol(:blasfeo_, flag, :trmv_ltu) + blasfeo_trmv_unn = Symbol(:blasfeo_, flag, :trmv_unn) + blasfeo_trmv_utn = Symbol(:blasfeo_, flag, :trmv_utn) + # blasfeo_trmv_unu = Symbol(:blasfeo_, flag, :trmv_unu) + # blasfeo_trmv_utu = Symbol(:blasfeo_, flag, :trmv_utu) + + @eval function Base.:*(A::LowerTriangular{$El, $Mat}, x::$Vec) + z = similar(x) + return mul!(z,A,x) + end + + @eval function LinearAlgebra.mul!(Y::$Vec, A::LowerTriangular{$El, $Mat}, B::$Vec) + @boundscheck begin + size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) + size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) + end + + $blasfeo_trmv_lnn( + size(A,1), + A.data, 0, 0, + B, 0, + Y, 0, + ) + return Y + end + + @eval function Base.:*(A::LowerTriangular{$El, Transpose{$El,$Mat}}, x::$Vec) + z = similar(x) + return mul!(z,A,x) + end + + @eval function LinearAlgebra.mul!(Y::$Vec, A::LowerTriangular{$El, Transpose{$El,$Mat}}, B::$Vec) + @boundscheck begin + size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) + size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) + end + $blasfeo_trmv_utn( + size(A,1), + A.data.parent, 0, 0, + B, 0, + Y, 0, + ) + return Y + end + + if El == :Cdouble # TODO(@anton) blasfeo_strmv_unn and blasfeo_strmv_lnu are unimplemented :( + @eval function Base.:*(A::UpperTriangular{$El, $Mat}, x::$Vec) + z = similar(x) + return mul!(z,A,x) + end + + @eval function LinearAlgebra.mul!(Y::$Vec, A::UpperTriangular{$El, $Mat}, B::$Vec) + @boundscheck begin + size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) + size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) + end + + $blasfeo_trmv_unn( + size(A,1), + A.data, 0, 0, + B, 0, + Y, 0, + ) + return Y + end + + @eval function Base.:*(A::UpperTriangular{$El, Transpose{$El,$Mat}}, x::$Vec) + z = similar(x) + return mul!(z,A,x) + end + + @eval function LinearAlgebra.mul!(Y::$Vec, A::UpperTriangular{$El, Transpose{$El,$Mat}}, B::$Vec) + @boundscheck begin + size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) + size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) + end + $blasfeo_trmv_ltn( + size(A,1), + A.data.parent, 0, 0, + B, 0, + Y, 0, + ) + return Y + end + + @eval function Base.:*(A::UnitLowerTriangular{$El, $Mat}, x::$Vec) + z = similar(x) + return mul!(z,A,x) + end + + @eval function LinearAlgebra.mul!(Y::$Vec, A::UnitLowerTriangular{$El, $Mat}, B::$Vec) + @boundscheck begin + size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) + size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) + end + + $blasfeo_trmv_lnu( + size(A,1), + A.data, 0, 0, + B, 0, + Y, 0, + ) + return Y + end + + @eval function Base.:*(A::UnitUpperTriangular{$El, Transpose{$El,$Mat}}, x::$Vec) + z = similar(x) + return mul!(z,A,x) + end + + @eval function LinearAlgebra.mul!(Y::$Vec, A::UnitUpperTriangular{$El, Transpose{$El,$Mat}}, B::$Vec) + @boundscheck begin + size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) + size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) + end + $blasfeo_trmv_ltu( + size(A,1), + A.data.parent, 0, 0, + B, 0, + Y, 0, + ) + return Y + end + end + + # need to implement trmv_unu + # @eval function Base.:*(A::UnitUpperTriangular{$El, $Mat}, x::$Vec) + # z = similar(x) + # return mul!(z,A,x) + # end + + # @eval function LinearAlgebra.mul!(Y::$Vec, A::UnitUpperTriangular{$El, $Mat}, B::$Vec) + # @boundscheck begin + # size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) + # size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) + # end + + # $blasfeo_trmv_unu( + # size(A,1), + # 1.0, A.data, 0, 0, + # B, 0, + # 0.0, B, 0, + # Y, 0, + # ) + # return Y + # end +end diff --git a/test/mat/level2.jl b/test/mat/level2.jl index cf62fe4..efb6dbd 100644 --- a/test/mat/level2.jl +++ b/test/mat/level2.jl @@ -1,84 +1,106 @@ @testset "Level 2 BLAS Operations" begin @testset for (MAT,VEC) in ((BlasfeoDmat,BlasfeoDvec), (BlasfeoSmat,BlasfeoSvec)) - A_orig = rand(eltype(VEC), 100, 100) - B_orig = rand(eltype(VEC), 95, 95) - C_orig = rand(eltype(VEC), 100, 50) - a_orig = rand(eltype(VEC), 100) - y_orig = rand(eltype(VEC), 100) - b_orig = rand(eltype(VEC), 90) - c_orig = rand(eltype(VEC), 50) - A = copy(A_orig) - B = copy(B_orig) - C = copy(C_orig) - a = copy(a_orig) - y = copy(y_orig) - b = copy(b_orig) - c = copy(c_orig) - A_blasfeo = MAT(A) - B_blasfeo = MAT(B) - C_blasfeo = MAT(C) - a_blasfeo = VEC(a) - y_blasfeo = VEC(y) - b_blasfeo = VEC(b) - c_blasfeo = VEC(c) + @testset for N in (10,11,91,100) + A_orig = rand(eltype(VEC), N, N) + B_orig = rand(eltype(VEC), N-5, N-5) + C_orig = rand(eltype(VEC), N, Int(round(N/2))) + a_orig = rand(eltype(VEC), N) + y_orig = rand(eltype(VEC), N) + b_orig = rand(eltype(VEC), N-10) + c_orig = rand(eltype(VEC), Int(round(N/2))) + A = copy(A_orig) + B = copy(B_orig) + C = copy(C_orig) + a = copy(a_orig) + y = copy(y_orig) + b = copy(b_orig) + c = copy(c_orig) + A_blasfeo = MAT(A) + B_blasfeo = MAT(B) + C_blasfeo = MAT(C) + a_blasfeo = VEC(a) + y_blasfeo = VEC(y) + b_blasfeo = VEC(b) + c_blasfeo = VEC(c) - @test_throws DimensionMismatch A_blasfeo*b_blasfeo - @test_throws DimensionMismatch B_blasfeo*a_blasfeo - @test_throws DimensionMismatch mul!(y_blasfeo, B_blasfeo,a_blasfeo) - @test_throws DimensionMismatch mul!(y_blasfeo, A_blasfeo,b_blasfeo) - @test_throws DimensionMismatch mul!(b_blasfeo, A_blasfeo,a_blasfeo) + @test_throws DimensionMismatch A_blasfeo*b_blasfeo + @test_throws DimensionMismatch B_blasfeo*a_blasfeo + @test_throws DimensionMismatch mul!(y_blasfeo, B_blasfeo,a_blasfeo) + @test_throws DimensionMismatch mul!(y_blasfeo, A_blasfeo,b_blasfeo) + @test_throws DimensionMismatch mul!(b_blasfeo, A_blasfeo,a_blasfeo) - # Test nondestructively - @test A*a ≈ A_blasfeo*a_blasfeo - @test isa(A_blasfeo*a_blasfeo, VEC) - @test transpose(A)*a ≈ transpose(A_blasfeo)*a_blasfeo - @test isa(transpose(A_blasfeo)*a_blasfeo, VEC) + # Test nondestructively + @test A*a ≈ A_blasfeo*a_blasfeo + @test isa(A_blasfeo*a_blasfeo, VEC) + @test transpose(A)*a ≈ transpose(A_blasfeo)*a_blasfeo + @test isa(transpose(A_blasfeo)*a_blasfeo, VEC) - # test mul! - @test mul!(y,A,a) ≈ mul!(y_blasfeo, A_blasfeo, a_blasfeo) - @test mul!(y,transpose(A),a) ≈ mul!(y_blasfeo, transpose(A_blasfeo), a_blasfeo) + # test mul! + @test mul!(y,A,a) ≈ mul!(y_blasfeo, A_blasfeo, a_blasfeo) + @test mul!(y,transpose(A),a) ≈ mul!(y_blasfeo, transpose(A_blasfeo), a_blasfeo) - # test asymmetric - @test C*c ≈ C_blasfeo*c_blasfeo - @test isa(C_blasfeo*c_blasfeo, VEC) - @test transpose(C)*y ≈ transpose(C_blasfeo)*y_blasfeo - @test isa(transpose(C_blasfeo)*y_blasfeo, VEC) + # test asymmetric + @test C*c ≈ C_blasfeo*c_blasfeo + @test isa(C_blasfeo*c_blasfeo, VEC) + @test transpose(C)*y ≈ transpose(C_blasfeo)*y_blasfeo + @test isa(transpose(C_blasfeo)*y_blasfeo, VEC) - # test symmetric - A_sym_L = Symmetric(A,:L) - A_sym_L_blasfeo = Symmetric(A_blasfeo,:L) - A_sym_U = Symmetric(A,:U) - A_sym_U_blasfeo = Symmetric(A_blasfeo,:U) + # test symmetric + A_sym_L = Symmetric(A,:L) + A_sym_L_blasfeo = Symmetric(A_blasfeo,:L) + A_sym_U = Symmetric(A,:U) + A_sym_U_blasfeo = Symmetric(A_blasfeo,:U) - @test A_sym_L*a ≈ A_sym_L_blasfeo*a_blasfeo - @test isa(A_sym_L_blasfeo*a_blasfeo, VEC) - @test A_sym_U*a ≈ A_sym_U_blasfeo*a_blasfeo - @test isa(A_sym_U_blasfeo*a_blasfeo, VEC) + @test A_sym_L*a ≈ A_sym_L_blasfeo*a_blasfeo + @test isa(A_sym_L_blasfeo*a_blasfeo, VEC) + @test A_sym_U*a ≈ A_sym_U_blasfeo*a_blasfeo + @test isa(A_sym_U_blasfeo*a_blasfeo, VEC) - # test triangular - A_lnn = LowerTriangular(A) - A_lnn_blasfeo = LowerTriangular(A_blasfeo) - A_unn = UpperTriangular(A) - A_unn_blasfeo = UpperTriangular(A_blasfeo) - A_lnu = UnitLowerTriangular(A) - A_lnu_blasfeo = UnitLowerTriangular(A_blasfeo) - A_unu = UnitUpperTriangular(A) - A_unu_blasfeo = UnitUpperTriangular(A_blasfeo) + # test triangular + A_lnn = LowerTriangular(A) + A_lnn_blasfeo = LowerTriangular(A_blasfeo) + A_unn = UpperTriangular(A) + A_unn_blasfeo = UpperTriangular(A_blasfeo) + A_lnu = UnitLowerTriangular(A) + A_lnu_blasfeo = UnitLowerTriangular(A_blasfeo) + A_unu = UnitUpperTriangular(A) + A_unu_blasfeo = UnitUpperTriangular(A_blasfeo) - # test lnn - @test A_lnn*a ≈ A_lnn_blasfeo*a_blasfeo - @test isa(A_lnn_blasfeo*a_blasfeo, VEC) - if MAT == BlasfeoDmat # TODO(@anton) not implemented upstream - # test unn - @test A_unn*a ≈ A_unn_blasfeo*a_blasfeo - @test isa(A_unn_blasfeo*a_blasfeo, VEC) - # test lnu - @test A_lnu*a ≈ A_lnu_blasfeo*a_blasfeo - @test isa(A_lnu_blasfeo*a_blasfeo, VEC) - # test unu - # TODO(@anton) not implemented upstream - #@test A_unu*a ≈ A_unu_blasfeo*a_blasfeo - #@test isa(A_unu_blasfeo*a_blasfeo, VEC) + # test lower triangular + @test A_lnn*a ≈ A_lnn_blasfeo*a_blasfeo + @test isa(A_lnn_blasfeo*a_blasfeo, VEC) + + # test upper triangular transpose + @test transpose(A_unn)*a ≈ transpose(A_unn_blasfeo)*a_blasfeo + @test isa(transpose(A_unn_blasfeo)*a_blasfeo, VEC) + + if MAT == BlasfeoDmat # TODO(@anton) some things not implemented upstream yet + # test lower triangular transpose + @test transpose(A_lnn)*a ≈ transpose(A_lnn_blasfeo)*a_blasfeo + @test isa(transpose(A_lnn_blasfeo)*a_blasfeo, VEC) + + # test upper triangular + @test A_unn*a ≈ A_unn_blasfeo*a_blasfeo + @test isa(A_unn_blasfeo*a_blasfeo, VEC) + + # test lower unit triangular + @test A_lnu*a ≈ A_lnu_blasfeo*a_blasfeo + @test isa(A_lnu_blasfeo*a_blasfeo, VEC) + + # test lower unit triangular transpose + @test transpose(A_lnu)*a ≈ transpose(A_lnu_blasfeo)*a_blasfeo + @test isa(transpose(A_lnu_blasfeo)*a_blasfeo, VEC) + + # test upper unit triangular + # TODO(@anton) not implemented upstream + #@test A_unu*a ≈ A_unu_blasfeo*a_blasfeo + #@test isa(A_unu_blasfeo*a_blasfeo, VEC) + + # test upper unit triangular transpose + # TODO(@anton) not implemented upstream + #@test transpose(A_unu)*a ≈ transpose(A_unu_blasfeo)*a_blasfeo + #@test isa(transpose(A_unu_blasfeo)*a_blasfeo, VEC) + end end end end From 6fba4b07e969e84972e545a9fbe9910d4f380df1 Mon Sep 17 00:00:00 2001 From: Anton Pozharskiy Date: Sun, 23 Aug 2026 15:46:11 +0200 Subject: [PATCH 5/9] only test implemented methods for single precision --- src/level2/trmv.jl | 35 +++++++++++++++++------------------ test/mat/level2.jl | 8 ++++---- 2 files changed, 21 insertions(+), 22 deletions(-) diff --git a/src/level2/trmv.jl b/src/level2/trmv.jl index 461dd89..c801aa6 100644 --- a/src/level2/trmv.jl +++ b/src/level2/trmv.jl @@ -32,26 +32,25 @@ for (El, Mat, Vec, flag) in [ return Y end - @eval function Base.:*(A::LowerTriangular{$El, Transpose{$El,$Mat}}, x::$Vec) - z = similar(x) - return mul!(z,A,x) - end - - @eval function LinearAlgebra.mul!(Y::$Vec, A::LowerTriangular{$El, Transpose{$El,$Mat}}, B::$Vec) - @boundscheck begin - size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) - size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) + if El == :Cdouble # TODO(@anton) blasfeo_strmv_unn and blasfeo_strmv_lnu are unimplemented :( + @eval function Base.:*(A::LowerTriangular{$El, Transpose{$El,$Mat}}, x::$Vec) + z = similar(x) + return mul!(z,A,x) end - $blasfeo_trmv_utn( - size(A,1), - A.data.parent, 0, 0, - B, 0, - Y, 0, - ) - return Y - end - if El == :Cdouble # TODO(@anton) blasfeo_strmv_unn and blasfeo_strmv_lnu are unimplemented :( + @eval function LinearAlgebra.mul!(Y::$Vec, A::LowerTriangular{$El, Transpose{$El,$Mat}}, B::$Vec) + @boundscheck begin + size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) + size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) + end + $blasfeo_trmv_utn( + size(A,1), + A.data.parent, 0, 0, + B, 0, + Y, 0, + ) + return Y + end @eval function Base.:*(A::UpperTriangular{$El, $Mat}, x::$Vec) z = similar(x) return mul!(z,A,x) diff --git a/test/mat/level2.jl b/test/mat/level2.jl index efb6dbd..c08634e 100644 --- a/test/mat/level2.jl +++ b/test/mat/level2.jl @@ -69,12 +69,12 @@ # test lower triangular @test A_lnn*a ≈ A_lnn_blasfeo*a_blasfeo @test isa(A_lnn_blasfeo*a_blasfeo, VEC) - - # test upper triangular transpose - @test transpose(A_unn)*a ≈ transpose(A_unn_blasfeo)*a_blasfeo - @test isa(transpose(A_unn_blasfeo)*a_blasfeo, VEC) if MAT == BlasfeoDmat # TODO(@anton) some things not implemented upstream yet + # test upper triangular transpose + @test transpose(A_unn)*a ≈ transpose(A_unn_blasfeo)*a_blasfeo + @test isa(transpose(A_unn_blasfeo)*a_blasfeo, VEC) + # test lower triangular transpose @test transpose(A_lnn)*a ≈ transpose(A_lnn_blasfeo)*a_blasfeo @test isa(transpose(A_lnn_blasfeo)*a_blasfeo, VEC) From 9f6090f8103378ce19efeb7a63b8049355213190 Mon Sep 17 00:00:00 2001 From: Anton Pozharskiy Date: Sun, 23 Aug 2026 19:25:43 +0200 Subject: [PATCH 6/9] fix indent --- src/level2/trmv.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/level2/trmv.jl b/src/level2/trmv.jl index c801aa6..2e65930 100644 --- a/src/level2/trmv.jl +++ b/src/level2/trmv.jl @@ -12,7 +12,7 @@ for (El, Mat, Vec, flag) in [ # blasfeo_trmv_unu = Symbol(:blasfeo_, flag, :trmv_unu) # blasfeo_trmv_utu = Symbol(:blasfeo_, flag, :trmv_utu) - @eval function Base.:*(A::LowerTriangular{$El, $Mat}, x::$Vec) + @eval function Base.:*(A::LowerTriangular{$El, $Mat}, x::$Vec) z = similar(x) return mul!(z,A,x) end From 3aa148f64017454f0af82aeecc45d47cfff3b330 Mon Sep 17 00:00:00 2001 From: Anton Pozharskiy Date: Sun, 23 Aug 2026 19:25:59 +0200 Subject: [PATCH 7/9] add diagonal specialization --- src/level2.jl | 1 + src/level2/dimv.jl | 27 +++++++++++++++++++++++++++ test/mat/level2.jl | 6 ++++++ 3 files changed, 34 insertions(+) create mode 100644 src/level2/dimv.jl diff --git a/src/level2.jl b/src/level2.jl index b824d08..ca371c7 100644 --- a/src/level2.jl +++ b/src/level2.jl @@ -1,3 +1,4 @@ include("level2/gemv.jl") include("level2/symv.jl") include("level2/trmv.jl") +include("level2/dimv.jl") diff --git a/src/level2/dimv.jl b/src/level2/dimv.jl new file mode 100644 index 0000000..5ee5095 --- /dev/null +++ b/src/level2/dimv.jl @@ -0,0 +1,27 @@ +for (El, Mat, Vec, flag) in [ + (:Cdouble, :BlasfeoDmat, :BlasfeoDvec, :d), + (:Cfloat, :BlasfeoSmat, :BlasfeoSvec, :s), + ] + + blasfeo_vecmul = Symbol(:blasfeo_, flag, :vecmul) + + @eval function Base.:*(A::Diagonal{$El, $Vec}, x::$Vec) + z = similar(x) + return mul!(z,A,x) + end + + @eval function LinearAlgebra.mul!(Y::$Vec, A::Diagonal{$El, $Vec}, B::$Vec) + @boundscheck begin + size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) + size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) + end + + $blasfeo_vecmul( + size(A,1), + A.diag, 0, + B, 0, + Y, 0, + ) + return Y + end +end diff --git a/test/mat/level2.jl b/test/mat/level2.jl index c08634e..c1f53b5 100644 --- a/test/mat/level2.jl +++ b/test/mat/level2.jl @@ -101,6 +101,12 @@ #@test transpose(A_unu)*a ≈ transpose(A_unu_blasfeo)*a_blasfeo #@test isa(transpose(A_unu_blasfeo)*a_blasfeo, VEC) end + + # test diagonal + Y_diag = Diagonal(y) + Y_diag_blasfeo = Diagonal(y_blasfeo) + @test Y_diag*a ≈ Y_diag_blasfeo*a_blasfeo + @test isa(Y_diag_blasfeo*a_blasfeo, VEC) end end end From 5d50d7ec5f6c3af3d53a5d3335df488efc0b0242 Mon Sep 17 00:00:00 2001 From: Anton Pozharskiy Date: Sun, 23 Aug 2026 21:33:00 +0200 Subject: [PATCH 8/9] add triangular solve for already supported solves --- src/level2.jl | 1 + src/level2/trmv.jl | 13 ++++ src/level2/trsv.jl | 143 ++++++++++++++++++++++++++++++++++++++++ test/mat/level2.jl | 159 +++++++++++++++++++++++++++++---------------- 4 files changed, 261 insertions(+), 55 deletions(-) create mode 100644 src/level2/trsv.jl diff --git a/src/level2.jl b/src/level2.jl index ca371c7..9cc1a3f 100644 --- a/src/level2.jl +++ b/src/level2.jl @@ -2,3 +2,4 @@ include("level2/gemv.jl") include("level2/symv.jl") include("level2/trmv.jl") include("level2/dimv.jl") +include("level2/trsv.jl") diff --git a/src/level2/trmv.jl b/src/level2/trmv.jl index 2e65930..dfb9217 100644 --- a/src/level2/trmv.jl +++ b/src/level2/trmv.jl @@ -12,6 +12,7 @@ for (El, Mat, Vec, flag) in [ # blasfeo_trmv_unu = Symbol(:blasfeo_, flag, :trmv_unu) # blasfeo_trmv_utu = Symbol(:blasfeo_, flag, :trmv_utu) + #-------------------------------------------------------------------------------------------------------------------------------# @eval function Base.:*(A::LowerTriangular{$El, $Mat}, x::$Vec) z = similar(x) return mul!(z,A,x) @@ -31,8 +32,10 @@ for (El, Mat, Vec, flag) in [ ) return Y end + #-------------------------------------------------------------------------------------------------------------------------------# if El == :Cdouble # TODO(@anton) blasfeo_strmv_unn and blasfeo_strmv_lnu are unimplemented :( + #-------------------------------------------------------------------------------------------------------------------------------# @eval function Base.:*(A::LowerTriangular{$El, Transpose{$El,$Mat}}, x::$Vec) z = similar(x) return mul!(z,A,x) @@ -51,6 +54,9 @@ for (El, Mat, Vec, flag) in [ ) return Y end + #-------------------------------------------------------------------------------------------------------------------------------# + + #-------------------------------------------------------------------------------------------------------------------------------# @eval function Base.:*(A::UpperTriangular{$El, $Mat}, x::$Vec) z = similar(x) return mul!(z,A,x) @@ -70,7 +76,9 @@ for (El, Mat, Vec, flag) in [ ) return Y end + #-------------------------------------------------------------------------------------------------------------------------------# + #-------------------------------------------------------------------------------------------------------------------------------# @eval function Base.:*(A::UpperTriangular{$El, Transpose{$El,$Mat}}, x::$Vec) z = similar(x) return mul!(z,A,x) @@ -89,7 +97,9 @@ for (El, Mat, Vec, flag) in [ ) return Y end + #-------------------------------------------------------------------------------------------------------------------------------# + #-------------------------------------------------------------------------------------------------------------------------------# @eval function Base.:*(A::UnitLowerTriangular{$El, $Mat}, x::$Vec) z = similar(x) return mul!(z,A,x) @@ -109,7 +119,9 @@ for (El, Mat, Vec, flag) in [ ) return Y end + #-------------------------------------------------------------------------------------------------------------------------------# + #-------------------------------------------------------------------------------------------------------------------------------# @eval function Base.:*(A::UnitUpperTriangular{$El, Transpose{$El,$Mat}}, x::$Vec) z = similar(x) return mul!(z,A,x) @@ -128,6 +140,7 @@ for (El, Mat, Vec, flag) in [ ) return Y end + #-------------------------------------------------------------------------------------------------------------------------------# end # need to implement trmv_unu diff --git a/src/level2/trsv.jl b/src/level2/trsv.jl new file mode 100644 index 0000000..0ab3ff5 --- /dev/null +++ b/src/level2/trsv.jl @@ -0,0 +1,143 @@ +for (El, Mat, Vec, flag) in [ + (:Cdouble, :BlasfeoDmat, :BlasfeoDvec, :d), + (:Cfloat, :BlasfeoSmat, :BlasfeoSvec, :s), + ] + + blasfeo_trsv_lnn = Symbol(:blasfeo_, flag, :trsv_lnn) + blasfeo_trsv_ltn = Symbol(:blasfeo_, flag, :trsv_ltn) + blasfeo_trsv_lnu = Symbol(:blasfeo_, flag, :trsv_lnu) + blasfeo_trsv_ltu = Symbol(:blasfeo_, flag, :trsv_ltu) + blasfeo_trsv_unn = Symbol(:blasfeo_, flag, :trsv_unn) + blasfeo_trsv_utn = Symbol(:blasfeo_, flag, :trsv_utn) + #blasfeo_trsv_unu = Symbol(:blasfeo_, flag, :trsv_unu) + #blasfeo_trsv_utu = Symbol(:blasfeo_, flag, :trsv_utu) + + #-------------------------------------------------------------------------------------------------------------------------------# + @eval function Base.:\(A::LowerTriangular{$El, $Mat}, x::$Vec) + z = similar(x) + return ldiv!(z,A,x) + end + + @eval function LinearAlgebra.ldiv!(Y::$Vec, A::LowerTriangular{$El, $Mat}, B::$Vec) + @boundscheck begin + size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) + size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) + end + + $blasfeo_trsv_lnn( + size(A,1), + A.data, 0, 0, + B, 0, + Y, 0, + ) + return Y + end + #-------------------------------------------------------------------------------------------------------------------------------# + + #-------------------------------------------------------------------------------------------------------------------------------# + @eval function Base.:\(A::LowerTriangular{$El, Transpose{$El,$Mat}}, x::$Vec) + z = similar(x) + return ldiv!(z,A,x) + end + + @eval function LinearAlgebra.ldiv!(Y::$Vec, A::LowerTriangular{$El, Transpose{$El,$Mat}}, B::$Vec) + @boundscheck begin + size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) + size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) + end + $blasfeo_trsv_utn( + size(A,1), + A.data.parent, 0, 0, + B, 0, + Y, 0, + ) + return Y + end + + #-------------------------------------------------------------------------------------------------------------------------------# + @eval function Base.:\(A::UpperTriangular{$El, Transpose{$El,$Mat}}, x::$Vec) + z = similar(x) + return ldiv!(z,A,x) + end + + @eval function LinearAlgebra.ldiv!(Y::$Vec, A::UpperTriangular{$El, Transpose{$El,$Mat}}, B::$Vec) + @boundscheck begin + size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) + size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) + end + $blasfeo_trsv_ltn( + size(A,1), + A.data.parent, 0, 0, + B, 0, + Y, 0, + ) + return Y + end + + if El == :Cdouble + + #-------------------------------------------------------------------------------------------------------------------------------# + @eval function Base.:\(A::UpperTriangular{$El, $Mat}, x::$Vec) + z = similar(x) + return ldiv!(z,A,x) + end + @eval function LinearAlgebra.ldiv!(Y::$Vec, A::UpperTriangular{$El, $Mat}, B::$Vec) + @boundscheck begin + size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) + size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) + end + + $blasfeo_trsv_unn( + size(A,1), + A.data, 0, 0, + B, 0, + Y, 0, + ) + return Y + end + #-------------------------------------------------------------------------------------------------------------------------------# + + #-------------------------------------------------------------------------------------------------------------------------------# + @eval function Base.:\(A::UnitLowerTriangular{$El, $Mat}, x::$Vec) + z = similar(x) + return ldiv!(z,A,x) + end + + @eval function LinearAlgebra.ldiv!(Y::$Vec, A::UnitLowerTriangular{$El, $Mat}, B::$Vec) + @boundscheck begin + size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) + size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) + end + + $blasfeo_trsv_lnu( + size(A,1), + A.data, 0, 0, + B, 0, + Y, 0, + ) + return Y + end + + #-------------------------------------------------------------------------------------------------------------------------------# + @eval function Base.:\(A::UnitUpperTriangular{$El, Transpose{$El,$Mat}}, x::$Vec) + z = similar(x) + return ldiv!(z,A,x) + end + + @eval function LinearAlgebra.ldiv!(Y::$Vec, A::UnitUpperTriangular{$El, Transpose{$El,$Mat}}, B::$Vec) + @boundscheck begin + size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) + size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) + end + $blasfeo_trsv_ltu( + size(A,1), + A.data.parent, 0, 0, + B, 0, + Y, 0, + ) + return Y + end + #-------------------------------------------------------------------------------------------------------------------------------# + end + # TODO(@anton) need to implement trsv_utu and trsv_unu upstream +end diff --git a/test/mat/level2.jl b/test/mat/level2.jl index c1f53b5..b60fddd 100644 --- a/test/mat/level2.jl +++ b/test/mat/level2.jl @@ -29,33 +29,39 @@ @test_throws DimensionMismatch mul!(y_blasfeo, A_blasfeo,b_blasfeo) @test_throws DimensionMismatch mul!(b_blasfeo, A_blasfeo,a_blasfeo) - # Test nondestructively - @test A*a ≈ A_blasfeo*a_blasfeo - @test isa(A_blasfeo*a_blasfeo, VEC) - @test transpose(A)*a ≈ transpose(A_blasfeo)*a_blasfeo - @test isa(transpose(A_blasfeo)*a_blasfeo, VEC) - - # test mul! - @test mul!(y,A,a) ≈ mul!(y_blasfeo, A_blasfeo, a_blasfeo) - @test mul!(y,transpose(A),a) ≈ mul!(y_blasfeo, transpose(A_blasfeo), a_blasfeo) - - # test asymmetric - @test C*c ≈ C_blasfeo*c_blasfeo - @test isa(C_blasfeo*c_blasfeo, VEC) - @test transpose(C)*y ≈ transpose(C_blasfeo)*y_blasfeo - @test isa(transpose(C_blasfeo)*y_blasfeo, VEC) - - # test symmetric + @testset "gemm" begin + # Test nondestructively + @test A*a ≈ A_blasfeo*a_blasfeo + @test isa(A_blasfeo*a_blasfeo, VEC) + @test transpose(A)*a ≈ transpose(A_blasfeo)*a_blasfeo + @test isa(transpose(A_blasfeo)*a_blasfeo, VEC) + + # test mul! + @test mul!(y,A,a) ≈ mul!(y_blasfeo, A_blasfeo, a_blasfeo) + @test mul!(y,transpose(A),a) ≈ mul!(y_blasfeo, transpose(A_blasfeo), a_blasfeo) + + # test nonsquare + @test C*c ≈ C_blasfeo*c_blasfeo + @test isa(C_blasfeo*c_blasfeo, VEC) + @test transpose(C)*y ≈ transpose(C_blasfeo)*y_blasfeo + @test isa(transpose(C_blasfeo)*y_blasfeo, VEC) + + end + + A_sym_L = Symmetric(A,:L) A_sym_L_blasfeo = Symmetric(A_blasfeo,:L) A_sym_U = Symmetric(A,:U) A_sym_U_blasfeo = Symmetric(A_blasfeo,:U) + @testset "symv" begin + # test symmetric - @test A_sym_L*a ≈ A_sym_L_blasfeo*a_blasfeo - @test isa(A_sym_L_blasfeo*a_blasfeo, VEC) - @test A_sym_U*a ≈ A_sym_U_blasfeo*a_blasfeo - @test isa(A_sym_U_blasfeo*a_blasfeo, VEC) + @test A_sym_L*a ≈ A_sym_L_blasfeo*a_blasfeo + @test isa(A_sym_L_blasfeo*a_blasfeo, VEC) + @test A_sym_U*a ≈ A_sym_U_blasfeo*a_blasfeo + @test isa(A_sym_U_blasfeo*a_blasfeo, VEC) + end # test triangular A_lnn = LowerTriangular(A) A_lnn_blasfeo = LowerTriangular(A_blasfeo) @@ -66,47 +72,90 @@ A_unu = UnitUpperTriangular(A) A_unu_blasfeo = UnitUpperTriangular(A_blasfeo) - # test lower triangular - @test A_lnn*a ≈ A_lnn_blasfeo*a_blasfeo - @test isa(A_lnn_blasfeo*a_blasfeo, VEC) - - if MAT == BlasfeoDmat # TODO(@anton) some things not implemented upstream yet - # test upper triangular transpose - @test transpose(A_unn)*a ≈ transpose(A_unn_blasfeo)*a_blasfeo - @test isa(transpose(A_unn_blasfeo)*a_blasfeo, VEC) + @testset "trmv" begin + # test lower triangular + @test A_lnn*a ≈ A_lnn_blasfeo*a_blasfeo + @test isa(A_lnn_blasfeo*a_blasfeo, VEC) - # test lower triangular transpose - @test transpose(A_lnn)*a ≈ transpose(A_lnn_blasfeo)*a_blasfeo - @test isa(transpose(A_lnn_blasfeo)*a_blasfeo, VEC) - - # test upper triangular - @test A_unn*a ≈ A_unn_blasfeo*a_blasfeo - @test isa(A_unn_blasfeo*a_blasfeo, VEC) + if MAT == BlasfeoDmat # TODO(@anton) some things not implemented upstream yet + # test upper triangular transpose + @test transpose(A_unn)*a ≈ transpose(A_unn_blasfeo)*a_blasfeo + @test isa(transpose(A_unn_blasfeo)*a_blasfeo, VEC) + + # test lower triangular transpose + @test transpose(A_lnn)*a ≈ transpose(A_lnn_blasfeo)*a_blasfeo + @test isa(transpose(A_lnn_blasfeo)*a_blasfeo, VEC) + + # test upper triangular + @test A_unn*a ≈ A_unn_blasfeo*a_blasfeo + @test isa(A_unn_blasfeo*a_blasfeo, VEC) + + # test lower unit triangular + @test A_lnu*a ≈ A_lnu_blasfeo*a_blasfeo + @test isa(A_lnu_blasfeo*a_blasfeo, VEC) + + # test lower unit triangular transpose + @test transpose(A_lnu)*a ≈ transpose(A_lnu_blasfeo)*a_blasfeo + @test isa(transpose(A_lnu_blasfeo)*a_blasfeo, VEC) + + # test upper unit triangular + # TODO(@anton) not implemented upstream + #@test A_unu*a ≈ A_unu_blasfeo*a_blasfeo + #@test isa(A_unu_blasfeo*a_blasfeo, VEC) + + # test upper unit triangular transpose + # TODO(@anton) not implemented upstream + #@test transpose(A_unu)*a ≈ transpose(A_unu_blasfeo)*a_blasfeo + #@test isa(transpose(A_unu_blasfeo)*a_blasfeo, VEC) + end + end + + @testset "trsv" begin + # test lower triangular + @test A_lnn\a ≈ A_lnn_blasfeo\a_blasfeo + @test isa(A_lnn_blasfeo\a_blasfeo, VEC) + - # test lower unit triangular - @test A_lnu*a ≈ A_lnu_blasfeo*a_blasfeo - @test isa(A_lnu_blasfeo*a_blasfeo, VEC) - - # test lower unit triangular transpose - @test transpose(A_lnu)*a ≈ transpose(A_lnu_blasfeo)*a_blasfeo - @test isa(transpose(A_lnu_blasfeo)*a_blasfeo, VEC) - - # test upper unit triangular - # TODO(@anton) not implemented upstream - #@test A_unu*a ≈ A_unu_blasfeo*a_blasfeo - #@test isa(A_unu_blasfeo*a_blasfeo, VEC) - - # test upper unit triangular transpose - # TODO(@anton) not implemented upstream - #@test transpose(A_unu)*a ≈ transpose(A_unu_blasfeo)*a_blasfeo - #@test isa(transpose(A_unu_blasfeo)*a_blasfeo, VEC) + if MAT == BlasfeoDmat # TODO(@anton) some things not implemented upstream yet + # test upper triangular + @test A_unn\a ≈ A_unn_blasfeo\a_blasfeo + @test isa(A_unn_blasfeo\a_blasfeo, VEC) + + # test lower triangular transpose + @test transpose(A_lnn)\a ≈ transpose(A_lnn_blasfeo)\a_blasfeo + @test isa(transpose(A_lnn_blasfeo)\a_blasfeo, VEC) + + # test upper triangular transpose + @test transpose(A_unn)\a ≈ transpose(A_unn_blasfeo)\a_blasfeo + @test isa(transpose(A_unn_blasfeo)\a_blasfeo, VEC) + + # test lower unit triangular + @test A_lnu\a ≈ A_lnu_blasfeo\a_blasfeo + @test isa(A_lnu_blasfeo\a_blasfeo, VEC) + + # test lower unit triangular transpose + @test transpose(A_lnu)\a ≈ transpose(A_lnu_blasfeo)\a_blasfeo + @test isa(transpose(A_lnu_blasfeo)\a_blasfeo, VEC) + + # test upper unit triangular + # TODO(@anton) not implemented upstream + #@test A_unu\a ≈ A_unu_blasfeo\a_blasfeo + #@test isa(A_unu_blasfeo\a_blasfeo, VEC) + + # test upper unit triangular transpose + # TODO(@anton) not implemented upstream + #@test transpose(A_unu)\a ≈ transpose(A_unu_blasfeo)\a_blasfeo + #@test isa(transpose(A_unu_blasfeo)\a_blasfeo, VEC) + end end # test diagonal Y_diag = Diagonal(y) Y_diag_blasfeo = Diagonal(y_blasfeo) - @test Y_diag*a ≈ Y_diag_blasfeo*a_blasfeo - @test isa(Y_diag_blasfeo*a_blasfeo, VEC) + @testset "dimv" begin + @test Y_diag*a ≈ Y_diag_blasfeo*a_blasfeo + @test isa(Y_diag_blasfeo*a_blasfeo, VEC) + end end end end From 8774537e30b4843ada85ba1a436c2179d06cae95 Mon Sep 17 00:00:00 2001 From: Anton Pozharskiy Date: Mon, 24 Aug 2026 08:35:48 +0200 Subject: [PATCH 9/9] add diagonal solve --- src/level2.jl | 3 ++- src/level2/disv.jl | 30 ++++++++++++++++++++++++++++++ test/mat/level2.jl | 4 ++++ 3 files changed, 36 insertions(+), 1 deletion(-) create mode 100644 src/level2/disv.jl diff --git a/src/level2.jl b/src/level2.jl index 9cc1a3f..423f6e4 100644 --- a/src/level2.jl +++ b/src/level2.jl @@ -1,5 +1,6 @@ include("level2/gemv.jl") include("level2/symv.jl") include("level2/trmv.jl") -include("level2/dimv.jl") include("level2/trsv.jl") +include("level2/dimv.jl") +include("level2/disv.jl") diff --git a/src/level2/disv.jl b/src/level2/disv.jl new file mode 100644 index 0000000..8ce688c --- /dev/null +++ b/src/level2/disv.jl @@ -0,0 +1,30 @@ +for (El, Mat, Vec, flag) in [ + (:Cdouble, :BlasfeoDmat, :BlasfeoDvec, :d), + (:Cfloat, :BlasfeoSmat, :BlasfeoSvec, :s), + ] + + blasfeo_vecmul = Symbol(:blasfeo_, flag, :vecmul) + + @eval function Base.:\(A::Diagonal{$El, $Vec}, x::$Vec) + z = similar(x) + return ldiv!(z,A,x) + end + + @eval function LinearAlgebra.ldiv!(Y::$Vec, A::Diagonal{$El, $Vec}, B::$Vec) + @boundscheck begin + size(A,2) == length(B) || throw(DimensionMismatch("Matrix second dimension doesn't match vector dimension")) + size(A,1) == length(Y) || throw(DimensionMismatch("Matrix first dimension doesn't match output vector dimension")) + end + + #TODO(@anton) make a blasfeo kernel for this! + Y .= $El(1.0)./A.diag + + $blasfeo_vecmul( + size(A,1), + B, 0, + Y, 0, + Y, 0, + ) + return Y + end +end diff --git a/test/mat/level2.jl b/test/mat/level2.jl index b60fddd..3f3856b 100644 --- a/test/mat/level2.jl +++ b/test/mat/level2.jl @@ -156,6 +156,10 @@ @test Y_diag*a ≈ Y_diag_blasfeo*a_blasfeo @test isa(Y_diag_blasfeo*a_blasfeo, VEC) end + @testset "disv" begin + @test Y_diag\a ≈ Y_diag_blasfeo\a_blasfeo + @test isa(Y_diag_blasfeo\a_blasfeo, VEC) + end end end end