From 58d5263f3d6f8177a14559f4779fcf72c1065b04 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Mon, 2 Sep 2024 08:47:59 +0530 Subject: [PATCH] Fix errant lmul! for tridiagonal and triangular (#55546) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit This method is incorrect. Firstly, the current implementation doesn't act in-place. Secondly, the result can't be stored in-place into a triangular matrix in general. This PR changes the implementation to convert the tridiagonal matrix to a `Diagonal` or a `Bidiagonal` one before attempting the `lmul!`. Currently, ```julia julia> T = Tridiagonal([1,2], [1,2,3], [1,2]); U = UpperTriangular(fill(2, 3, 3)); julia> lmul!(T, U) 3×3 Matrix{Int64}: 2 4 4 2 6 10 0 4 10 julia> U # not updated 3×3 UpperTriangular{Int64, Matrix{Int64}}: 2 2 2 ⋅ 2 2 ⋅ ⋅ 2 julia> parent(U) # except for the underlying storage 3×3 Matrix{Int64}: 2 2 2 0 2 2 0 0 2 ``` After this, ```julia julia> lmul!(T, U) ERROR: ArgumentError: matrix cannot be represented as Bidiagonal ``` I'm unsure if we want this method at all, since there isn't a corresponding `rmul!`, but I've left it there to avoid breakages, if any. --- stdlib/LinearAlgebra/src/triangular.jl | 2 -- stdlib/LinearAlgebra/test/triangular.jl | 2 -- 2 files changed, 4 deletions(-) diff --git a/stdlib/LinearAlgebra/src/triangular.jl b/stdlib/LinearAlgebra/src/triangular.jl index 923e13e488c85..8d32dac824ce8 100644 --- a/stdlib/LinearAlgebra/src/triangular.jl +++ b/stdlib/LinearAlgebra/src/triangular.jl @@ -953,8 +953,6 @@ isunit_char(::UnitUpperTriangular) = 'U' isunit_char(::LowerTriangular) = 'N' isunit_char(::UnitLowerTriangular) = 'U' -lmul!(A::Tridiagonal, B::AbstractTriangular) = A*full!(B) - # generic fallback for AbstractTriangular matrices outside of the four subtypes provided here _trimul!(C::AbstractVecOrMat, A::AbstractTriangular, B::AbstractVector) = lmul!(A, copyto!(C, B)) diff --git a/stdlib/LinearAlgebra/test/triangular.jl b/stdlib/LinearAlgebra/test/triangular.jl index 5f0a829f9cdda..b88be00b0209c 100644 --- a/stdlib/LinearAlgebra/test/triangular.jl +++ b/stdlib/LinearAlgebra/test/triangular.jl @@ -442,8 +442,6 @@ Base.getindex(A::MyTriangular, i::Int, j::Int) = A.data[i,j] debug && println("elty1: $elty1, A1: $t1, B: $eltyB") - Tri = Tridiagonal(rand(eltyB,n-1),rand(eltyB,n),rand(eltyB,n-1)) - @test lmul!(Tri,copy(A1)) ≈ Tri*M1 Tri = Tridiagonal(rand(eltyB,n-1),rand(eltyB,n),rand(eltyB,n-1)) C = Matrix{promote_type(elty1,eltyB)}(undef, n, n) mul!(C, Tri, A1)