improve copyto!(::AbstractMatrix,::SpecialMatrix)`

stdlib/LinearAlgebra/src/bidiag.jl

@@ -165,29 +165,34 @@ function Base.replace_in_print_matrix(A::Bidiagonal,i::Integer,j::Integer,s::Abs
 end
 
 #Converting from Bidiagonal to dense Matrix
-function Matrix{T}(A::Bidiagonal) where T
- n = size(A, 1)
- B = zeros(T, n, n)
- if n == 0
- return B
- end
- for i = 1:n - 1
- B[i,i] = A.dv[i]
- if A.uplo == 'U'
- B[i, i + 1] = A.ev[i]
- else
- B[i + 1, i] = A.ev[i]
- end
- end
- B[n,n] = A.dv[n]
- return B
-end
+Matrix{T}(A::Bidiagonal) where {T} = copyto!(Matrix{T}(undef, size(A)), A)
 Matrix(A::Bidiagonal{T}) where {T} = Matrix{T}(A)
 Array(A::Bidiagonal) = Matrix(A)
 promote_rule(::Type{Matrix{T}}, ::Type{<:Bidiagonal{S}}) where {T,S} =
  @isdefined(T) && @isdefined(S) ? Matrix{promote_type(T,S)} : Matrix
 promote_rule(::Type{Matrix}, ::Type{<:Bidiagonal}) = Matrix
 
+function copyto!(A::AbstractMatrix{T}, B::Bidiagonal) where {T}
+ require_one_based_indexing(A)
+ n = size(B, 1)
+ n == 0 && return A
+ if size(A) == (n, n)
+ fill!(A, zero(T))
+ @inbounds for i in 1:n - 1
+ A[i,i] = B.dv[i]
+ if B.uplo == 'U'
+ A[i, i + 1] = B.ev[i]
+ else
+ A[i + 1, i] = B.ev[i]
+ end
+ end
+ A[n,n] = B.dv[n]
+ return A
+ else
+ return @invoke copyto!(A::AbstractMatrix, B::AbstractMatrix)
+ end
+end
+
 #Converting from Bidiagonal to Tridiagonal
 function Tridiagonal{T}(A::Bidiagonal) where T
  dv = convert(AbstractVector{T}, A.dv)

stdlib/LinearAlgebra/src/diagonal.jl

@@ -91,6 +91,20 @@ similar(D::Diagonal, ::Type{T}) where {T} = Diagonal(similar(D.diag, T))
 similar(::Diagonal, ::Type{T}, dims::Union{Dims{1},Dims{2}}) where {T} = zeros(T, dims...)
 
 copyto!(D1::Diagonal, D2::Diagonal) = (copyto!(D1.diag, D2.diag); D1)
+function copyto!(A::AbstractMatrix{T}, D::Diagonal) where {T}
+ require_one_based_indexing(A)
+ n = length(D.diag)
+ n == 0 && return A
+ if size(A) == (n, n)
+ fill!(A, zero(T))
+ @inbounds for i in 1:n
+ A[i, i] = D.diag[i]
+ end
+ return A
+ else
+ return @invoke copyto!(A::AbstractMatrix, D::AbstractMatrix)
+ end
+end
 
 size(D::Diagonal) = (n = length(D.diag); (n,n))
 

stdlib/LinearAlgebra/src/special.jl

@@ -300,9 +300,9 @@ lmul!(Q::Adjoint{<:Any,<:QRPackedQ}, B::AbstractTriangular) = lmul!(Q, full!(B))
 function _qlmul(Q::AbstractQ, B)
  TQB = promote_type(eltype(Q), eltype(B))
  if size(Q.factors, 1) == size(B, 1)
- Bnew = Matrix{TQB}(B)
+ Bnew = copy_similar(B, TQB)
  elseif size(Q.factors, 2) == size(B, 1)
- Bnew = [Matrix{TQB}(B); zeros(TQB, size(Q.factors, 1) - size(B,1), size(B, 2))]
+ Bnew = [copy_similar(B, TQB); zeros(TQB, size(Q.factors, 1) - size(B,1), size(B, 2))]
  else
  throw(DimensionMismatch("first dimension of matrix must have size either $(size(Q.factors, 1)) or $(size(Q.factors, 2))"))
  end
@@ -331,9 +331,9 @@ function _qrmul(A, adjQ::Adjoint{<:Any,<:AbstractQ})
  Q = adjQ.parent
  TAQ = promote_type(eltype(A), eltype(Q))
  if size(A,2) == size(Q.factors, 1)
- Anew = Matrix{TAQ}(A)
+ Anew = copy_similar(A, TAQ)
  elseif size(A,2) == size(Q.factors,2)
- Anew = [Matrix{TAQ}(A) zeros(TAQ, size(A, 1), size(Q.factors, 1) - size(Q.factors, 2))]
+ Anew = [copy_similar(A, TAQ) zeros(TAQ, size(A, 1), size(Q.factors, 1) - size(Q.factors, 2))]
  else
  throw(DimensionMismatch("matrix A has dimensions $(size(A)) but matrix B has dimensions $(size(Q))"))
  end

stdlib/LinearAlgebra/src/tridiag.jl

@@ -122,17 +122,23 @@ SymTridiagonal(S::SymTridiagonal) = S
 AbstractMatrix{T}(S::SymTridiagonal) where {T} =
  SymTridiagonal(convert(AbstractVector{T}, S.dv)::AbstractVector{T},
  convert(AbstractVector{T}, S.ev)::AbstractVector{T})
-function Matrix{T}(M::SymTridiagonal) where T
+Matrix{T}(M::SymTridiagonal) where {T} = copyto!(Matrix{T}(undef, size(M)), M)
+function copyto!(A::AbstractMatrix{T}, M::SymTridiagonal) where {T}
+ require_one_based_indexing(A)
  n = size(M, 1)
- Mf = zeros(T, n, n)
- n == 0 && return Mf
- @inbounds for i = 1:n-1
- Mf[i,i] = symmetric(M.dv[i], :U)
- Mf[i+1,i] = transpose(M.ev[i])
- Mf[i,i+1] = M.ev[i]
+ n == 0 && return A
+ if size(A) == (n, n)
+ fill!(A, zero(T))
+ @inbounds for i in 1:n-1
+ A[i,i] = symmetric(M.dv[i], :U)
+ A[i+1,i] = transpose(M.ev[i])
+ A[i,i+1] = M.ev[i]
+ end
+ A[n,n] = symmetric(M.dv[n], :U)
+ return A
+ else
+ return @invoke copyto!(A::AbstractMatrix, M::AbstractMatrix)
  end
- Mf[n,n] = symmetric(M.dv[n], :U)
- return Mf
 end
 Matrix(M::SymTridiagonal{T}) where {T} = Matrix{T}(M)
 Array(M::SymTridiagonal) = Matrix(M)
@@ -571,16 +577,24 @@ function size(M::Tridiagonal, d::Integer)
  end
 end
 
-function Matrix{T}(M::Tridiagonal{T}) where T
- A = zeros(T, size(M))
- for i = 1:length(M.d)
- A[i,i] = M.d[i]
- end
- for i = 1:length(M.d)-1
- A[i+1,i] = M.dl[i]
- A[i,i+1] = M.du[i]
+Matrix{T}(M::Tridiagonal) where {T} = copyto!(Matrix{T}(undef, size(M)), M)
+
+function copyto!(A::AbstractMatrix{T}, M::Tridiagonal) where {T}
+ require_one_based_indexing(A)
+ n = size(M, 1)
+ n == 0 && return A
+ if size(A) == (n, n)
+ fill!(A, zero(T))
+ @inbounds for i = 1:n-1
+ A[i,i] = M.d[i]
+ A[i+1,i] = M.dl[i]
+ A[i,i+1] = M.du[i]
+ end
+ A[n,n] = M.d[n]
+ return A
+ else
+ @invoke copyto!(A::AbstractMatrix, M::AbstractMatrix)
  end
- A
 end
 Matrix(M::Tridiagonal{T}) where {T} = Matrix{T}(M)
 Array(M::Tridiagonal) = Matrix(M)

stdlib/LinearAlgebra/test/special.jl

@@ -429,12 +429,13 @@ end
  dl = [1, 1, 1]
  d = [1, 1, 1, 1]
  D = Diagonal(d)
- Bi = Bidiagonal(d, dl, :L)
+ Bl = Bidiagonal(d, dl, :L)
+ Bu = Bidiagonal(d, dl, :U)
  Tri = Tridiagonal(dl, d, dl)
  Sym = SymTridiagonal(d, dl)
  F = qr(ones(4, 1))
  A = F.Q'
- for A in (F.Q, F.Q'), B in (D, Bi, Tri, Sym)
+ for A in (F.Q, F.Q'), B in (D, Bl, Bu, Tri, Sym)
  @test B*A ≈ Matrix(B)*A
  @test A*B ≈ A*Matrix(B)
  end