Fix generic triangular solves with empty matrices (JuliaLang#54201)

Liozou · Apr 23, 2024 · 8945914 · 8945914
1 parent ef168a6
commit 8945914
Show file tree

Hide file tree

Showing 2 changed files with 15 additions and 0 deletions.
diff --git a/stdlib/LinearAlgebra/src/triangular.jl b/stdlib/LinearAlgebra/src/triangular.jl
@@ -1296,6 +1296,7 @@ function generic_trimatdiv!(C::AbstractVecOrMat, uploc, isunitc, tfun::Function,
  if size(C) != size(B)
  throw(DimensionMismatch(lazy"size of output, $(size(C)), does not match size of right hand side, $(size(B))"))
  end
+ iszero(mA) && return C
  oA = oneunit(eltype(A))
  @inbounds if uploc == 'U'
  if isunitc == 'N'
@@ -1432,6 +1433,7 @@ function generic_trimatdiv!(C::AbstractVecOrMat, uploc, isunitc, ::Function, xA:
  if size(C) != size(B)
  throw(DimensionMismatch(lazy"size of output, $(size(C)), does not match size of right hand side, $(size(B))"))
  end
+ iszero(mA) && return C
  oA = oneunit(eltype(A))
  @inbounds if uploc == 'U'
  if isunitc == 'N'

diff --git a/stdlib/LinearAlgebra/test/lu.jl b/stdlib/LinearAlgebra/test/lu.jl
@@ -486,4 +486,17 @@ end
  LinearAlgebra.generic_lufact!(fill(Inf, 2, 2), check=false)
 end
 
+@testset "lu for empty matrices" begin
+ for T in (Float64, BigFloat)
+ A = fill(T(0.0), 0, 0)
+ v = fill(T(1.0), 0, 10)
+ @test A \ v ≈ lu(A) \ v
+ vt = permutedims(v)
+ @test vt / A ≈ vt / lu(A)
+ B = UpperTriangular(transpose(fill(complex(T(0.0)), 0, 0)'))
+ @test B \ v ≈ v
+ @test vt / B ≈ vt
+ end
+end
+
 end # module TestLU