From 965b88440fbb74e03aeef29ad5d09111a48888eb Mon Sep 17 00:00:00 2001 From: Sacha Verweij Date: Mon, 16 Oct 2017 10:39:44 -0700 Subject: [PATCH] Make similar(S<:[Sym]Tridiagonal, [T,] shape...) return a sparse rather 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). --- base/linalg/lu.jl | 6 +++--- base/linalg/triangular.jl | 10 ++-------- base/linalg/tridiag.jl | 15 +++++++++++---- test/linalg/tridiag.jl | 6 ++++-- 4 files changed, 20 insertions(+), 17 deletions(-) diff --git a/base/linalg/lu.jl b/base/linalg/lu.jl index 086314a0c8eeee..9bae05f033af9b 100644 --- a/base/linalg/lu.jl +++ b/base/linalg/lu.jl @@ -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 diff --git a/base/linalg/triangular.jl b/base/linalg/triangular.jl index f0ab6b86ce9b4e..0157b7f8825be3 100644 --- a/base/linalg/triangular.jl +++ b/base/linalg/triangular.jl @@ -347,9 +347,7 @@ 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::UnitLowerTriangular) = ones(eltype(A), size(A,1)) -diag(A::UpperTriangular) = diag(A.data) diag(A::UnitUpperTriangular) = ones(eltype(A), size(A,1)) # Unary operations @@ -1430,13 +1428,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 diff --git a/base/linalg/tridiag.jl b/base/linalg/tridiag.jl index d84c53225b8e4d..bd992ac8dadadf 100644 --- a/base/linalg/tridiag.jl +++ b/base/linalg/tridiag.jl @@ -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)) @@ -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) diff --git a/test/linalg/tridiag.jl b/test/linalg/tridiag.jl index 8bda2ede296979..3d44af2bb437c6 100644 --- a/test/linalg/tridiag.jl +++ b/test/linalg/tridiag.jl @@ -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) @@ -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))