do not call ldlt! in the poly algorithm for hermitian matrices (#325)

Author: KristofferC

src/solvers/cholmod.jl

@@ -1577,12 +1577,7 @@ function \(A::RealHermSymComplexHermF64SSL, B::StridedVecOrMatInclAdjAndTrans)
     if issuccess(F)
         return \(F, B)
     else
-        ldlt!(F, A; check = false)
-        if issuccess(F)
-            return \(F, B)
-        else
-            return \(lu(SparseMatrixCSC{eltype(A), SuiteSparse_long}(A)), B)
-        end
+        return \(lu(SparseMatrixCSC{eltype(A), SuiteSparse_long}(A)), B)
     end
 end
 

test/cholmod.jl

@@ -927,6 +927,30 @@ end
     end
 end
 
+@testset "sym indefinite poly alg" begin
+    K = open(joinpath(@__DIR__, "matrices", "stiffness_sym_indef")) do io
+        ml = readline(io)
+        m = parse(Int, split(ml, "m = ")[2])
+        nl = readline(io)
+        n = parse(Int, split(nl, "n = ")[2])
+
+        colptrl = readline(io)
+        rowvall = readline(io)
+        nzvall = readline(io)
+
+        colptr = parse.(Int,     split(strip(split(colptrl, "colptr = ")[2], [']', '[']), ','))
+        rowval = parse.(Int,     split(strip(split(rowvall, "rowval = ")[2], [']', '[']), ','))
+        nzval =  parse.(Float64, split(strip(split(nzvall, "nzval = ")[2], [']', '[']), ','))
+
+        SparseMatrixCSC(m, n, colptr, rowval, nzval)
+    end
+
+    f = ones(size(K, 1))
+    u = K \ f
+    residual = norm(f - K * u) / norm(f)
+    @test residual < 1e-6 
+end
+
 end # Base.USE_GPL_LIBS
 
 end # module