Fix \ and .'\ for Symmetric and Hermitian sparse matrices. (#19242)

(cherry picked from commit d9d0fad)

base/linalg/factorization.jl

@@ -30,8 +30,7 @@ function (\){T<:BlasReal}(F::Factorization{T}, B::VecOrMat{Complex{T}})
 end
 
 for (f1, f2) in ((:\, :A_ldiv_B!),
-                 (:Ac_ldiv_B, :Ac_ldiv_B!),
-                 (:At_ldiv_B, :At_ldiv_B!))
+                 (:Ac_ldiv_B, :Ac_ldiv_B!))
     @eval begin
         function $f1(F::Factorization, B::AbstractVecOrMat)
             TFB = typeof(one(eltype(F)) / one(eltype(B)))
@@ -48,6 +47,10 @@ for f in (:A_ldiv_B!, :Ac_ldiv_B!, :At_ldiv_B!)
         $f(A, copy!(Y, B))
 end
 
+# fallback methods for transposed solves
+At_ldiv_B{T<:Real}(F::Factorization{T}, B::AbstractVecOrMat) = Ac_ldiv_B(F, B)
+At_ldiv_B(F::Factorization, B) = conj.(Ac_ldiv_B(F, conj.(B)))
+
 """
     A_ldiv_B!([Y,] A, B) -> Y
 

base/sparse/cholmod.jl

@@ -1498,6 +1498,19 @@ Ac_ldiv_B(L::Factor, B::VecOrMat) = convert(Matrix, solve(CHOLMOD_A, L, Dense(B)
 Ac_ldiv_B(L::Factor, B::Sparse) = spsolve(CHOLMOD_A, L, B)
 Ac_ldiv_B(L::Factor, B::SparseVecOrMat) = Ac_ldiv_B(L, Sparse(B))
 
+for f in (:\, :Ac_ldiv_B)
+    @eval function ($f)(A::Union{Symmetric{Float64,SparseMatrixCSC{Float64,SuiteSparse_long}},
+                          Hermitian{Float64,SparseMatrixCSC{Float64,SuiteSparse_long}},
+                          Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},SuiteSparse_long}}}, B::StridedVecOrMat)
+        try
+            return ($f)(cholfact(A), B)
+        catch e
+            isa(e, LinAlg.PosDefException) || rethrow(e)
+            return ($f)(ldltfact(A) , B)
+        end
+    end
+end
+
 ## Other convenience methods
 function diag{Tv}(F::Factor{Tv})
     f = unsafe_load(get(F.p))

base/sparse/linalg.jl

@@ -864,12 +864,7 @@ function (\)(A::SparseMatrixCSC, B::AbstractVecOrMat)
             return UpperTriangular(A) \ B
         end
         if ishermitian(A)
-            try
-                return cholfact(Hermitian(A)) \ B
-            catch e
-                isa(e, PosDefException) || rethrow(e)
-                return ldltfact(Hermitian(A)) \ B
-            end
+            Hermitian(A) \ B
         end
         return lufact(A) \ B
     else

test/linalg/lu.jl

@@ -44,7 +44,7 @@ for eltya in (Float32, Float64, Complex64, Complex128, BigFloat, Int)
         εb = eps(abs(float(one(eltyb))))
         ε = max(εa,εb)
 
-debug && println("(Automatic) Square LU decomposition")
+debug && println("(Automatic) Square LU decomposition. eltya: $eltya, eltyb: $eltyb")
         κ     = cond(a,1)
         lua   = factorize(a)
         @test_throws KeyError lua[:Z]

test/sparsedir/cholmod.jl

@@ -382,6 +382,7 @@ for elty in (Float64, Complex{Float64})
     @test_throws DimensionMismatch F\CHOLMOD.Sparse(sparse(ones(elty, 4)))
     @test F'\ones(elty, 5) ≈ full(A1pd)'\ones(5)
     @test F'\sparse(ones(elty, 5)) ≈ full(A1pd)'\ones(5)
+    @test F.'\ones(elty, 5) ≈ conj(A1pd)'\ones(elty, 5)
     @test logdet(F) ≈ logdet(full(A1pd))
     @test det(F) == exp(logdet(F))
     let # to test supernodal, we must use a larger matrix
@@ -671,3 +672,19 @@ let m = 400, n = 500
     @test s.is_super == 0
     @test F\b ≈ ones(m + n)
 end
+
+# Test that \ and '\ and .'\ work for Symmetric and Hermitian. This is just
+# a dispatch exercise so it doesn't matter that the complex matrix has
+# zero imaginary parts
+let Apre = sprandn(10, 10, 0.2) - I
+    for A in (Symmetric(Apre), Hermitian(Apre),
+              Symmetric(Apre + 10I), Hermitian(Apre + 10I),
+              Hermitian(complex(Apre)), Hermitian(complex(Apre) + 10I))
+
+        x = ones(10)
+        b = A*x
+        @test x ≈ A\b
+        @test A.'\b ≈ A'\b
+    end
+end
+