Skip to content

Commit

Permalink
Fix \ and .'\ for Symmetric and Hermitian sparse matrices.
Browse files Browse the repository at this point in the history
Undeprecate conj for sparse matrices.
  • Loading branch information
andreasnoack committed Nov 14, 2016
1 parent 735554d commit d433e6f
Show file tree
Hide file tree
Showing 6 changed files with 38 additions and 10 deletions.
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

0 comments on commit d433e6f

Please sign in to comment.