Make ldlt throw with zero pivot (#38051)

JuliaLang · Oct 29, 2020 · f613b55 · f613b55
1 parent 72971c4
commit f613b55
Show file tree

Hide file tree

Showing 3 changed files with 15 additions and 7 deletions.
diff --git a/stdlib/LinearAlgebra/src/dense.jl b/stdlib/LinearAlgebra/src/dense.jl
@@ -1284,12 +1284,10 @@ function factorize(A::StridedMatrix{T}) where T
                 return Bidiagonal(diag(A), diag(A, 1), :U)
             end
             if utri1
-                if (herm & (T <: Complex)) | sym
-                    try
-                        return ldlt!(SymTridiagonal(diag(A), diag(A, -1)))
-                    catch
-                    end
-                end
+                # TODO: enable once a specialized, non-dense bunchkaufman method exists
+                # if (herm & (T <: Complex)) | sym
+                    # return bunchkaufman(SymTridiagonal(diag(A), diag(A, -1)))
+                # end
                 return lu(Tridiagonal(diag(A, -1), diag(A), diag(A, 1)))
             end
         end

diff --git a/stdlib/LinearAlgebra/src/ldlt.jl b/stdlib/LinearAlgebra/src/ldlt.jl
@@ -115,7 +115,8 @@ function ldlt!(S::SymTridiagonal{T,V}) where {T,V}
     n = size(S,1)
     d = S.dv
     e = S.ev
-    @inbounds @simd for i = 1:n-1
+    @inbounds for i in 1:n-1
+        iszero(d[i]) && throw(ZeroPivotException(i))
         e[i] /= d[i]
         d[i+1] -= e[i]^2*d[i]
     end

diff --git a/stdlib/LinearAlgebra/test/lu.jl b/stdlib/LinearAlgebra/test/lu.jl
@@ -210,6 +210,15 @@ dimg  = randn(n)/2
             @test luhs.L*luhs.U ≈ luhs.P*Matrix(HS)
         end
     end
+
+    @testset "Factorization of symtridiagonal dense matrix with zero ldlt-pivot (#38026)" begin
+        A = [0.0 -1.0 0.0 0.0
+            -1.0 0.0 0.0 0.0
+            0.0 0.0 0.0 -1.0
+            0.0 0.0 -1.0 0.0]
+        F = factorize(A)
+        @test all((!isnan).(Matrix(F)))
+    end
 end
 
 @testset "Singular matrices" for T in (Float64, ComplexF64)