Merge pull request #24170 from Sacha0/similarstruc

make similar(A, [T,] shape...) for structured A yield a sparse rather than dense array

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 #

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)

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

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"))
 
 

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

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)

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

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

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))