Skip to content

Commit

Permalink
Make similar(S<:[Sym]Tridiagonal, [T,] shape...) return a sparse rath…
Browse files Browse the repository at this point in the history
…er than dense array.

Also fix a minor bug where lufact(A::Tridiagonal{T}) for !(T<:AbstractFloat) would dispatch
to lufact!(B, ...) for non-tridiagonal rather than tridiagonal B downstream.

Also make certain that diag(A::Union{UpperTriangular,LowerTriangular}, k) yields diag(A.data, k)
independent of k (i.e. particularly fro diag(A) = diag(A, k) as well).
  • Loading branch information
Sacha0 committed Oct 19, 2017
1 parent 41c97c9 commit 5962bc6
Show file tree
Hide file tree
Showing 4 changed files with 22 additions and 17 deletions.
6 changes: 3 additions & 3 deletions base/linalg/lu.jl
Original file line number Diff line number Diff line change
Expand Up @@ -140,20 +140,20 @@ true
"""
function lufact(A::AbstractMatrix{T}, pivot::Union{Val{false}, Val{true}}) where T
S = typeof(zero(T)/one(T))
AA = similar(A, S, size(A))
AA = similar(A, S)
copy!(AA, A)
lufact!(AA, pivot)
end
# We can't assume an ordered field so we first try without pivoting
function lufact(A::AbstractMatrix{T}) where T
S = typeof(zero(T)/one(T))
AA = similar(A, S, size(A))
AA = similar(A, S)
copy!(AA, A)
F = lufact!(AA, Val(false))
if issuccess(F)
return F
else
AA = similar(A, S, size(A))
AA = similar(A, S)
copy!(AA, A)
return lufact!(AA, Val(true))
end
Expand Down
12 changes: 4 additions & 8 deletions base/linalg/triangular.jl
Original file line number Diff line number Diff line change
Expand Up @@ -347,9 +347,9 @@ adjoint!(A::UnitLowerTriangular) = UnitUpperTriangular(copytri!(A.data, 'L' , tr
adjoint!(A::UpperTriangular) = LowerTriangular(copytri!(A.data, 'U' , true))
adjoint!(A::UnitUpperTriangular) = UnitLowerTriangular(copytri!(A.data, 'U' , true))

diag(A::LowerTriangular) = diag(A.data)
diag(A::LowerTriangular, k::Integer = 0) = diag(A.data, k)
diag(A::UnitLowerTriangular) = ones(eltype(A), size(A,1))
diag(A::UpperTriangular) = diag(A.data)
diag(A::UpperTriangular, k::Integer = 0) = diag(A.data, k)
diag(A::UnitUpperTriangular) = ones(eltype(A), size(A,1))

# Unary operations
Expand Down Expand Up @@ -1430,13 +1430,9 @@ end
## the element type doesn't have to be stable under division whereas that is
## necessary in the general triangular solve problem.

## Some Triangular-Triangular cases. We might want to write taylored methods
## Some Triangular-Triangular cases. We might want to write tailored methods
## for these cases, but I'm not sure it is worth it.
for t in (UpperTriangular, UnitUpperTriangular, LowerTriangular, UnitLowerTriangular)
@eval begin
(*)(A::Tridiagonal, B::$t) = A_mul_B!(Matrix(A), B)
end
end
(*)(A::Union{Tridiagonal,SymTridiagonal}, B::AbstractTriangular) = A_mul_B!(Matrix(A), B)

for (f1, f2) in ((:*, :A_mul_B!), (:\, :A_ldiv_B!))
@eval begin
Expand Down
15 changes: 11 additions & 4 deletions base/linalg/tridiag.jl
Original file line number Diff line number Diff line change
Expand Up @@ -108,7 +108,11 @@ function size(A::SymTridiagonal, d::Integer)
end
end

similar(S::SymTridiagonal, ::Type{T}) where {T} = SymTridiagonal{T}(similar(S.dv, T), similar(S.ev, T))
# For S<:SymTridiagonal, similar(S[, neweltype]) should yield a SymTridiagonal matrix.
# On the other hand, similar(S, [neweltype,] shape...) should yield a sparse matrix.
# The first method below effects the former, and the second the latter.
similar(S::SymTridiagonal, ::Type{T}) where {T} = SymTridiagonal(similar(S.dv, T), similar(S.ev, T))
similar(S::SymTridiagonal, ::Type{T}, dims::Union{Dims{1},Dims{2}}) where {T} = spzeros(T, dims...)

#Elementary operations
broadcast(::typeof(abs), M::SymTridiagonal) = SymTridiagonal(abs.(M.dv), abs.(M.ev))
Expand Down Expand Up @@ -499,9 +503,12 @@ end
convert(::Type{Matrix}, M::Tridiagonal{T}) where {T} = convert(Matrix{T}, M)
convert(::Type{Array}, M::Tridiagonal) = convert(Matrix, M)
full(M::Tridiagonal) = convert(Array, M)
function similar(M::Tridiagonal, ::Type{T}) where T
Tridiagonal{T}(similar(M.dl, T), similar(M.d, T), similar(M.du, T))
end

# For M<:Tridiagonal, similar(M[, neweltype]) should yield a Tridiagonal matrix.
# On the other hand, similar(M, [neweltype,] shape...) should yield a sparse matrix.
# The first method below effects the former, and the second the latter.
similar(M::Tridiagonal, ::Type{T}) where {T} = Tridiagonal(similar(M.dl, T), similar(M.d, T), similar(M.du, T))
similar(M::Tridiagonal, ::Type{T}, dims::Union{Dims{1},Dims{2}}) where {T} = spzeros(T, dims...)

# Operations on Tridiagonal matrices
copy!(dest::Tridiagonal, src::Tridiagonal) = (copy!(dest.dl, src.dl); copy!(dest.d, src.d); copy!(dest.du, src.du); dest)
Expand Down
6 changes: 4 additions & 2 deletions test/linalg/tridiag.jl
Original file line number Diff line number Diff line change
Expand Up @@ -124,7 +124,8 @@ guardsrand(123) do
end
@test isa(similar(A), mat_type{elty})
@test isa(similar(A, Int), mat_type{Int})
@test isa(similar(A, Int, (3, 2)), Matrix{Int})
@test isa(similar(A, (3, 2)), SparseMatrixCSC)
@test isa(similar(A, Int, (3, 2)), SparseMatrixCSC{Int})
@test size(A, 3) == 1
@test size(A, 1) == n
@test size(A) == (n, n)
Expand Down Expand Up @@ -289,7 +290,8 @@ guardsrand(123) do
@testset "similar" begin
@test isa(similar(Ts), SymTridiagonal{elty})
@test isa(similar(Ts, Int), SymTridiagonal{Int})
@test isa(similar(Ts, Int, (3,2)), Matrix{Int})
@test isa(similar(Ts, (3, 2)), SparseMatrixCSC)
@test isa(similar(Ts, Int, (3, 2)), SparseMatrixCSC{Int})
end

@test first(logabsdet(Tldlt)) first(logabsdet(Fs))
Expand Down

0 comments on commit 5962bc6

Please sign in to comment.