Merge pull request #15198 from JuliaLang/anj/similar

Use Matrix for similar(SpecialMatrix, shape) instead of throwing
JuliaLang · Feb 24, 2016 · 4d8f5bf · 4d8f5bf
2 parents 86f42bc + e696f0a
commit 4d8f5bf
Show file tree

Hide file tree

Showing 11 changed files with 56 additions and 28 deletions.
diff --git a/base/abstractarray.jl b/base/abstractarray.jl
@@ -189,7 +189,7 @@ end
 ## Constructors ##
 
 # default arguments to similar()
-similar{T}(a::AbstractArray{T})                          = similar(a, T, size(a))
+similar{T}(a::AbstractArray{T})                          = similar(a, T)
 similar(   a::AbstractArray, T::Type)                    = similar(a, T, size(a))
 similar{T}(a::AbstractArray{T}, dims::DimsInteger)       = similar(a, T, dims)
 similar{T}(a::AbstractArray{T}, dims::Integer...)        = similar(a, T, dims)

diff --git a/base/linalg/bidiag.jl b/base/linalg/bidiag.jl
@@ -5,7 +5,7 @@ type Bidiagonal{T} <: AbstractMatrix{T}
     dv::Vector{T} # diagonal
     ev::Vector{T} # sub/super diagonal
     isupper::Bool # is upper bidiagonal (true) or lower (false)
-    function Bidiagonal{T}(dv::Vector{T}, ev::Vector{T}, isupper::Bool)
+    function Bidiagonal(dv::Vector{T}, ev::Vector{T}, isupper::Bool)
         if length(ev) != length(dv)-1
             throw(DimensionMismatch("length of diagonal vector is $(length(dv)), length of off-diagonal vector is $(length(ev))"))
         end
@@ -101,6 +101,8 @@ convert{Tnew,Told}(::Type{AbstractMatrix{Tnew}}, A::Bidiagonal{Told}) = convert(
 
 big(B::Bidiagonal) = Bidiagonal(big(B.dv), big(B.ev), B.isupper)
 
+similar{T}(B::Bidiagonal, ::Type{T}) = Bidiagonal{T}(similar(B.dv, T), similar(B.ev, T), B.isupper)
+
 ###################
 # LAPACK routines #
 ###################

diff --git a/base/linalg/diagonal.jl b/base/linalg/diagonal.jl
@@ -14,11 +14,8 @@ convert{T}(::Type{AbstractMatrix{T}}, D::Diagonal) = convert(Diagonal{T}, D)
 convert{T}(::Type{UpperTriangular}, A::Diagonal{T}) = UpperTriangular(A)
 convert{T}(::Type{LowerTriangular}, A::Diagonal{T}) = LowerTriangular(A)
 
-function similar{T}(D::Diagonal, ::Type{T}, d::Tuple{Int,Int})
-    if d[1] != d[2]
-        throw(ArgumentError("diagonal matrix must be square"))
-    end
-    return Diagonal{T}(Array(T,d[1]))
+function similar{T}(D::Diagonal, ::Type{T})
+    return Diagonal{T}(similar(D.diag, T))
 end
 
 copy!(D1::Diagonal, D2::Diagonal) = (copy!(D1.diag, D2.diag); D1)

diff --git a/base/linalg/symmetric.jl b/base/linalg/symmetric.jl
@@ -33,6 +33,18 @@ end
     @inbounds r = (A.uplo == 'U') == (i < j) ? A.data[i, j] : conj(A.data[j, i])
     r
 end
+
+similar{T}(A::Symmetric, ::Type{T}) = Symmetric(similar(A.data, T))
+# Hermitian version can be simplified when check for imaginary part of
+# diagonal in Hermitian has been removed
+function similar{T}(A::Hermitian, ::Type{T})
+    B = similar(A.data, T)
+    for i = 1:size(A,1)
+        B[i,i] = 0
+    end
+    return Hermitian(B)
+end
+
 full(A::Symmetric) = copytri!(copy(A.data), A.uplo)
 full(A::Hermitian) = copytri!(copy(A.data), A.uplo, true)
 parent(A::HermOrSym) = A.data

diff --git a/base/linalg/triangular.jl b/base/linalg/triangular.jl
@@ -24,14 +24,8 @@ for t in (:LowerTriangular, :UnitLowerTriangular, :UpperTriangular,
         convert{Tnew,Told,S}(::Type{AbstractMatrix{Tnew}}, A::$t{Told,S}) = convert($t{Tnew}, A)
         convert{T,S}(::Type{Matrix}, A::$t{T,S}) = convert(Matrix{T}, A)
 
-        function similar{T,S,Tnew}(A::$t{T,S}, ::Type{Tnew}, dims::Dims)
-            if length(dims) != 2
-                throw(ArgumentError("Triangular matrix must have two dimensions"))
-            end
-            if dims[1] != dims[2]
-                throw(ArgumentError("Triangular matrix must be square"))
-            end
-            B = similar(A.data, Tnew, dims)
+        function similar{T,S,Tnew}(A::$t{T,S}, ::Type{Tnew})
+            B = similar(A.data, Tnew)
             return $t(B)
         end
 

diff --git a/base/linalg/tridiag.jl b/base/linalg/tridiag.jl
@@ -60,6 +60,8 @@ function size(A::SymTridiagonal, d::Integer)
     end
 end
 
+similar{T}(S::SymTridiagonal, ::Type{T}) = SymTridiagonal{T}(similar(S.dv, T), similar(S.ev, T))
+
 #Elementary operations
 for func in (:conj, :copy, :round, :trunc, :floor, :ceil, :abs, :real, :imag)
     @eval ($func)(M::SymTridiagonal) = SymTridiagonal(($func)(M.dv), ($func)(M.ev))
@@ -361,11 +363,8 @@ function convert{T}(::Type{Matrix{T}}, M::Tridiagonal{T})
     A
 end
 convert{T}(::Type{Matrix}, M::Tridiagonal{T}) = convert(Matrix{T}, M)
-function similar(M::Tridiagonal, T, dims::Dims)
-    if length(dims) != 2 || dims[1] != dims[2]
-        throw(DimensionMismatch("Tridiagonal matrices must be square"))
-    end
-    Tridiagonal{T}(similar(M.dl), similar(M.d), similar(M.du), similar(M.du2))
+function similar{T}(M::Tridiagonal, ::Type{T})
+    Tridiagonal{T}(similar(M.dl, T), similar(M.d, T), similar(M.du, T), similar(M.du2, T))
 end
 
 # Operations on Tridiagonal matrices

diff --git a/test/linalg/bidiag.jl b/test/linalg/bidiag.jl
@@ -32,7 +32,7 @@ for relty in (Int, Float32, Float64, BigFloat), elty in (relty, Complex{relty})
     @test_throws DimensionMismatch Bidiagonal(dv,ones(elty,n),true)
     @test_throws ArgumentError Bidiagonal(dv,ev)
 
-    debug && println("getindex, setindex!, and size")
+    debug && println("getindex, setindex!, size, and similar")
     BD = Bidiagonal(dv, ev, true)
     @test_throws BoundsError BD[n+1,1]
     @test BD[2,2] == dv[2]
@@ -44,6 +44,9 @@ for relty in (Int, Float32, Float64, BigFloat), elty in (relty, Complex{relty})
     @test BD == cBD
     @test_throws ArgumentError size(BD,0)
     @test size(BD,3) == 1
+    @test isa(similar(BD), Bidiagonal{elty})
+    @test isa(similar(BD, Int), Bidiagonal{Int})
+    @test isa(similar(BD, Int, (3,2)), Matrix{Int})
 
     debug && println("show")
     dstring = sprint(Base.print_matrix,BD.dv')

diff --git a/test/linalg/diagonal.jl b/test/linalg/diagonal.jl
@@ -158,8 +158,10 @@ for relty in (Float32, Float64, BigFloat), elty in (relty, Complex{relty})
             end
 
             #similar
-            @test_throws ArgumentError similar(D, eltype(D), (n,n+1))
-            @test length(diag(similar(D, eltype(D), (n,n)))) == n
+            @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})
 
             #10036
             @test issym(D2)

diff --git a/test/linalg/symmetric.jl b/test/linalg/symmetric.jl
@@ -46,6 +46,16 @@ let n=10
 
         debug && println("\ntype of a: ", eltya, "\n")
 
+        # similar
+        @test isa(similar(Symmetric(asym)), Symmetric{eltya})
+        @test isa(similar(Hermitian(asym)), Hermitian{eltya})
+        @test isa(similar(Symmetric(asym), Int), Symmetric{Int})
+        @test isa(similar(Hermitian(asym), Int), Hermitian{Int})
+        @test isa(similar(Symmetric(asym), (3,2)), Matrix{eltya})
+        @test isa(similar(Hermitian(asym), (3,2)), Matrix{eltya})
+        @test isa(similar(Symmetric(asym), Int, (3,2)), Matrix{Int})
+        @test isa(similar(Hermitian(asym), Int, (3,2)), Matrix{Int})
+
         # full
         @test asym == full(Hermitian(asym))
 

diff --git a/test/linalg/triangular.jl b/test/linalg/triangular.jl
@@ -38,14 +38,18 @@ for elty1 in (Float32, Float64, BigFloat, Complex64, Complex128, Complex{BigFloa
 
         # similar
         @test isa(similar(A1), t1)
-        @test_throws ArgumentError similar(A1,typeof(A1),(n,n+1))
-        @test_throws ArgumentError similar(A1,typeof(A1),(n,n,n))
+        @test eltype(similar(A1)) == elty1
+        @test isa(similar(A1, Int), t1)
+        @test eltype(similar(A1, Int)) == Int
+        @test isa(similar(A1, (3,2)), Matrix{elty1})
+        @test isa(similar(A1, Int, (3,2)), Matrix{Int})
 
         # getindex
         ## Linear indexing
         for i = 1:length(A1)
             @test A1[i] == full(A1)[i]
         end
+        @test isa(A1[2:4,1], Vector)
 
         ## Cartesian indexing
         for i = 1:size(A1, 1)

diff --git a/test/linalg/tridiag.jl b/test/linalg/tridiag.jl
@@ -121,6 +121,11 @@ for elty in (Float32, Float64, Complex64, Complex128, Int)
             @test_approx_eq x invFsv
             @test_approx_eq full(full(Tldlt)) Fs
         end
+
+        # similar
+        @test isa(similar(Ts), SymTridiagonal{elty})
+        @test isa(similar(Ts, Int), SymTridiagonal{Int})
+        @test isa(similar(Ts, Int, (3,2)), Matrix{Int})
     end
 
     # eigenvalues/eigenvectors of symmetric tridiagonal
@@ -365,9 +370,9 @@ let n = 12 #Size of matrix problem to test
         @test size(B) == size(A)
         copy!(B,A)
         @test B == A
-        @test_throws DimensionMismatch similar(A,(n,n,2))
-        @test_throws DimensionMismatch similar(A,(n+1,n))
-        @test_throws DimensionMismatch similar(A,(n,n+1))
+        @test isa(similar(A), Tridiagonal{elty})
+        @test isa(similar(A, Int), Tridiagonal{Int})
+        @test isa(similar(A, Int, (3,2)), Matrix{Int})
         @test size(A,3) == 1
         @test_throws ArgumentError size(A,0)