Skip to content

Commit

Permalink
Fix errant lmul! for tridiagonal and triangular (JuliaLang#55546)
Browse files Browse the repository at this point in the history
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.
  • Loading branch information
jishnub authored Sep 2, 2024
1 parent 4b99e99 commit 58d5263
Show file tree
Hide file tree
Showing 2 changed files with 0 additions and 4 deletions.
2 changes: 0 additions & 2 deletions stdlib/LinearAlgebra/src/triangular.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand Down
2 changes: 0 additions & 2 deletions stdlib/LinearAlgebra/test/triangular.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down

0 comments on commit 58d5263

Please sign in to comment.