diff --git a/base/linalg/bidiag.jl b/base/linalg/bidiag.jl index 46241dc6ef1fd..43f112a0e9411 100644 --- a/base/linalg/bidiag.jl +++ b/base/linalg/bidiag.jl @@ -170,7 +170,12 @@ convert(::Type{AbstractMatrix{T}}, A::Bidiagonal) where {T} = convert(Bidiagonal broadcast(::typeof(big), B::Bidiagonal) = Bidiagonal(big.(B.dv), big.(B.ev), B.uplo) -similar(B::Bidiagonal, ::Type{T}) where {T} = Bidiagonal{T}(similar(B.dv, T), similar(B.ev, T), B.uplo) +# For B<:Bidiagonal, similar(B[, neweltype]) should yield a Bidiagonal matrix. +# On the other hand, similar(B, [neweltype,] shape...) should yield a sparse matrix. +# The first method below effects the former, and the second the latter. +similar(B::Bidiagonal, ::Type{T}) where {T} = Bidiagonal(similar(B.dv, T), similar(B.ev, T), B.uplo) +similar(B::Bidiagonal, ::Type{T}, dims::Union{Dims{1},Dims{2}}) where {T} = spzeros(T, dims...) + ################### # LAPACK routines # diff --git a/base/linalg/diagonal.jl b/base/linalg/diagonal.jl index b8f0b53ebd7de..d0273215f8ffe 100644 --- a/base/linalg/diagonal.jl +++ b/base/linalg/diagonal.jl @@ -56,9 +56,11 @@ convert(::Type{Matrix}, D::Diagonal) = diagm(0 => D.diag) convert(::Type{Array}, D::Diagonal) = convert(Matrix, D) full(D::Diagonal) = convert(Array, D) -function similar(D::Diagonal, ::Type{T}) where T - return Diagonal{T}(similar(D.diag, T)) -end +# For D<:Diagonal, similar(D[, neweltype]) should yield a Diagonal matrix. +# On the other hand, similar(D, [neweltype,] shape...) should yield a sparse matrix. +# The first method below effects the former, and the second the latter. +similar(D::Diagonal, ::Type{T}) where {T} = Diagonal(similar(D.diag, T)) +similar(D::Diagonal, ::Type{T}, dims::Union{Dims{1},Dims{2}}) where {T} = spzeros(T, dims...) copy!(D1::Diagonal, D2::Diagonal) = (copy!(D1.diag, D2.diag); D1) @@ -305,6 +307,11 @@ end A_rdiv_Bc!(A::AbstractMatrix{T}, D::Diagonal{T}) where {T} = A_rdiv_B!(A, conj(D)) A_rdiv_Bt!(A::AbstractMatrix{T}, D::Diagonal{T}) where {T} = A_rdiv_B!(A, D) +(\)(F::Factorization, D::Diagonal) = + A_ldiv_B!(F, Matrix{typeof(oneunit(eltype(D))/oneunit(eltype(F)))}(D)) +Ac_ldiv_B(F::Factorization, D::Diagonal) = + Ac_ldiv_B!(F, Matrix{typeof(oneunit(eltype(D))/oneunit(eltype(F)))}(D)) + # Methods to resolve ambiguities with `Diagonal` @inline *(rowvec::RowVector, D::Diagonal) = transpose(D * transpose(rowvec)) @inline A_mul_Bt(D::Diagonal, rowvec::RowVector) = D*transpose(rowvec) diff --git a/base/linalg/lu.jl b/base/linalg/lu.jl index 086314a0c8eee..9bae05f033af9 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/special.jl b/base/linalg/special.jl index d3a4ab3b4482e..48a8053653d3d 100644 --- a/base/linalg/special.jl +++ b/base/linalg/special.jl @@ -52,13 +52,13 @@ convert(::Type{Diagonal}, A::AbstractTriangular) = isdiag(A) ? Diagonal(diag(A)) : throw(ArgumentError("matrix cannot be represented as Diagonal")) convert(::Type{Bidiagonal}, A::AbstractTriangular) = - isbanded(A, 0, 1) ? Bidiagonal(diag(A), diag(A, 1), :U) : # is upper bidiagonal - isbanded(A, -1, 0) ? Bidiagonal(diag(A), diag(A, -1), :L) : # is lower bidiagonal + isbanded(A, 0, 1) ? Bidiagonal(diag(A, 0), diag(A, 1), :U) : # is upper bidiagonal + isbanded(A, -1, 0) ? Bidiagonal(diag(A, 0), diag(A, -1), :L) : # is lower bidiagonal throw(ArgumentError("matrix cannot be represented as Bidiagonal")) convert(::Type{SymTridiagonal}, A::AbstractTriangular) = convert(SymTridiagonal, convert(Tridiagonal, A)) convert(::Type{Tridiagonal}, A::AbstractTriangular) = - isbanded(A, -1, 1) ? Tridiagonal(diag(A, -1), diag(A), diag(A, 1)) : # is tridiagonal + isbanded(A, -1, 1) ? Tridiagonal(diag(A, -1), diag(A, 0), diag(A, 1)) : # is tridiagonal throw(ArgumentError("matrix cannot be represented as Tridiagonal")) diff --git a/base/linalg/triangular.jl b/base/linalg/triangular.jl index 0be40f7afe69d..fa7d325d5d3c7 100644 --- a/base/linalg/triangular.jl +++ b/base/linalg/triangular.jl @@ -1434,11 +1434,7 @@ end ## Some Triangular-Triangular cases. We might want to write taylored 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 b1ebae7609966..f6f71b0f71b40 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/bidiag.jl b/test/linalg/bidiag.jl index 111dfc4f242de..6594ad02e40af 100644 --- a/test/linalg/bidiag.jl +++ b/test/linalg/bidiag.jl @@ -78,8 +78,11 @@ srand(1) @test size(ubd, 3) == 1 # bidiagonal similar @test isa(similar(ubd), Bidiagonal{elty}) + @test similar(ubd).uplo == ubd.uplo @test isa(similar(ubd, Int), Bidiagonal{Int}) - @test isa(similar(ubd, Int, (3, 2)), Matrix{Int}) + @test similar(ubd, Int).uplo == ubd.uplo + @test isa(similar(ubd, (3, 2)), SparseMatrixCSC) + @test isa(similar(ubd, Int, (3, 2)), SparseMatrixCSC{Int}) end @testset "show" begin diff --git a/test/linalg/diagonal.jl b/test/linalg/diagonal.jl index d2afb5280501f..0d6e12b694901 100644 --- a/test/linalg/diagonal.jl +++ b/test/linalg/diagonal.jl @@ -228,8 +228,8 @@ srand(1) @testset "similar" begin @test isa(similar(D), Diagonal{elty}) @test isa(similar(D, Int), Diagonal{Int}) - @test isa(similar(D, (3,2)), Matrix{elty}) - @test isa(similar(D, Int, (3,2)), Matrix{Int}) + @test isa(similar(D, (3,2)), SparseMatrixCSC{elty}) + @test isa(similar(D, Int, (3,2)), SparseMatrixCSC{Int}) end # Issue number 10036 diff --git a/test/linalg/tridiag.jl b/test/linalg/tridiag.jl index 437e67b9f877a..3649ab415e5b8 100644 --- a/test/linalg/tridiag.jl +++ b/test/linalg/tridiag.jl @@ -125,7 +125,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) @@ -290,7 +291,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))