Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Merged
merged 1 commit into from
Nov 16, 2016
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion base/deprecated.jl
Original file line number Diff line number Diff line change
Expand Up @@ -855,7 +855,7 @@ end
for f in (:sin, :sinh, :sind, :asin, :asinh, :asind,
:tan, :tanh, :tand, :atan, :atanh, :atand,
:sinpi, :cosc, :ceil, :floor, :trunc, :round, :real, :imag,
:log1p, :expm1, :abs, :abs2, :conj,
:log1p, :expm1, :abs, :abs2,
:log, :log2, :log10, :exp, :exp2, :exp10, :sinc, :cospi,
:cos, :cosh, :cosd, :acos, :acosd,
:cot, :coth, :cotd, :acot, :acotd,
Expand Down
7 changes: 5 additions & 2 deletions base/linalg/factorization.jl
Original file line number Diff line number Diff line change
Expand Up @@ -39,8 +39,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)))
Expand All @@ -57,6 +56,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

Expand Down
13 changes: 13 additions & 0 deletions base/sparse/cholmod.jl
Original file line number Diff line number Diff line change
Expand Up @@ -1504,6 +1504,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))
Expand Down
7 changes: 1 addition & 6 deletions base/sparse/linalg.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion test/linalg/lu.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down
17 changes: 17 additions & 0 deletions test/sparse/cholmod.jl
Original file line number Diff line number Diff line change
Expand Up @@ -381,6 +381,7 @@ for elty in (Float64, Complex{Float64})
@test_throws DimensionMismatch F\CHOLMOD.Sparse(sparse(ones(elty, 4)))
@test F'\ones(elty, 5) ≈ Array(A1pd)'\ones(5)
@test F'\sparse(ones(elty, 5)) ≈ Array(A1pd)'\ones(5)
@test F.'\ones(elty, 5) ≈ conj(A1pd)'\ones(elty, 5)
@test logdet(F) ≈ logdet(Array(A1pd))
@test det(F) == exp(logdet(F))
let # to test supernodal, we must use a larger matrix
Expand Down Expand Up @@ -670,3 +671,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