Use copyto! in converting Diagonal/Bidiagonal/Tridiagonal to …

…`Matrix` (#53912)

With this, we may convert a structured matrix to `Matrix` even if its
`eltype` doesn't support `zero(T)`, as long as we may index into the
matrix and the elements have `zero` defined for themselves. This makes
the following work:
```julia
julia> D = Diagonal(fill(Diagonal([1,3]), 2))
2×2 Diagonal{Diagonal{Int64, Vector{Int64}}, Vector{Diagonal{Int64, Vector{Int64}}}}:
 [1 0; 0 3]      ⋅     
     ⋅       [1 0; 0 3]

julia> Matrix{eltype(D)}(D)
2×2 Matrix{Diagonal{Int64, Vector{Int64}}}:
 [1 0; 0 3]  [0 0; 0 0]
 [0 0; 0 0]  [1 0; 0 3]
```
We also may materialize partly initialized matrices:
```julia
julia> D = Diagonal(Vector{BigInt}(undef, 2))
2×2 Diagonal{BigInt, Vector{BigInt}}:
 #undef    ⋅
   ⋅     #undef

julia> Matrix{eltype(D)}(D)
2×2 Matrix{BigInt}:
 #undef    0
   0     #undef
```

The performance seems identical for numeric matrices.

stdlib/LinearAlgebra/src/bidiag.jl

@@ -195,19 +195,14 @@ end
 
 #Converting from Bidiagonal to dense Matrix
 function Matrix{T}(A::Bidiagonal) where T
-    n = size(A, 1)
-    B = Matrix{T}(undef, n, n)
-    n == 0 && return B
-    n > 1 && fill!(B, zero(T))
-    @inbounds 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
+    B = Matrix{T}(undef, size(A))
+    if haszero(T) # optimized path for types with zero(T) defined
+        size(B,1) > 1 && fill!(B, zero(T))
+        copyto!(view(B, diagind(B)), A.dv)
+        copyto!(view(B, diagind(B, A.uplo == 'U' ? 1 : -1)), A.ev)
+    else
+        copyto!(B, A)
     end
-    B[n,n] = A.dv[n]
     return B
 end
 Matrix(A::Bidiagonal{T}) where {T} = Matrix{promote_type(T, typeof(zero(T)))}(A)

stdlib/LinearAlgebra/src/diagonal.jl

@@ -116,11 +116,12 @@ AbstractMatrix{T}(D::Diagonal{T}) where {T} = copy(D)
 Matrix(D::Diagonal{T}) where {T} = Matrix{promote_type(T, typeof(zero(T)))}(D)
 Array(D::Diagonal{T}) where {T} = Matrix(D)
 function Matrix{T}(D::Diagonal) where {T}
-    n = size(D, 1)
-    B = Matrix{T}(undef, n, n)
-    n > 1 && fill!(B, zero(T))
-    @inbounds for i in 1:n
-        B[i,i] = D.diag[i]
+    B = Matrix{T}(undef, size(D))
+    if haszero(T) # optimized path for types with zero(T) defined
+        size(B,1) > 1 && fill!(B, zero(T))
+        copyto!(view(B, diagind(B)), D.diag)
+    else
+        copyto!(B, D)
     end
     return B
 end

stdlib/LinearAlgebra/src/tridiag.jl

@@ -135,13 +135,17 @@ function Matrix{T}(M::SymTridiagonal) where T
     n = size(M, 1)
     Mf = Matrix{T}(undef, n, n)
     n == 0 && return Mf
-    n > 2 && fill!(Mf, zero(T))
-    @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]
+    if haszero(T) # optimized path for types with zero(T) defined
+        n > 2 && fill!(Mf, zero(T))
+        @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]
+        end
+        Mf[n,n] = symmetric(M.dv[n], :U)
+    else
+        copyto!(Mf, M)
     end
-    Mf[n,n] = symmetric(M.dv[n], :U)
     return Mf
 end
 Matrix(M::SymTridiagonal{T}) where {T} = Matrix{promote_type(T, typeof(zero(T)))}(M)
@@ -586,15 +590,14 @@ axes(M::Tridiagonal) = (ax = axes(M.d,1); (ax, ax))
 
 function Matrix{T}(M::Tridiagonal) where {T}
     A = Matrix{T}(undef, size(M))
-    n = length(M.d)
-    n == 0 && return A
-    n > 2 && fill!(A, zero(T))
-    for i in 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]
+    if haszero(T) # optimized path for types with zero(T) defined
+        size(A,1) > 2 && fill!(A, zero(T))
+        copyto!(view(A, diagind(A)), M.d)
+        copyto!(view(A, diagind(A,1)), M.du)
+        copyto!(view(A, diagind(A,-1)), M.dl)
+    else
+        copyto!(A, M)
+    end
     A
 end
 Matrix(M::Tridiagonal{T}) where {T} = Matrix{promote_type(T, typeof(zero(T)))}(M)

stdlib/LinearAlgebra/test/bidiag.jl

@@ -884,4 +884,17 @@ end
     @test mul!(C1, B, sv, 1, 2) == mul!(C2, B, v, 1 ,2)
 end
 
+@testset "Matrix conversion for non-numeric and undef" begin
+    B = Bidiagonal(Vector{BigInt}(undef, 4), fill(big(3), 3), :U)
+    M = Matrix(B)
+    B[diagind(B)] .= 4
+    M[diagind(M)] .= 4
+    @test diag(B) == diag(M)
+
+    B = Bidiagonal(fill(Diagonal([1,3]), 3), fill(Diagonal([1,3]), 2), :U)
+    M = Matrix{eltype(B)}(B)
+    @test M isa Matrix{eltype(B)}
+    @test M == B
+end
+
 end # module TestBidiagonal

stdlib/LinearAlgebra/test/diagonal.jl

@@ -1292,4 +1292,17 @@ end
     @test yadj == x'
 end
 
+@testset "Matrix conversion for non-numeric and undef" begin
+    D = Diagonal(Vector{BigInt}(undef, 4))
+    M = Matrix(D)
+    D[diagind(D)] .= 4
+    M[diagind(M)] .= 4
+    @test diag(D) == diag(M)
+
+    D = Diagonal(fill(Diagonal([1,3]), 2))
+    M = Matrix{eltype(D)}(D)
+    @test M isa Matrix{eltype(D)}
+    @test M == D
+end
+
 end # module TestDiagonal

stdlib/LinearAlgebra/test/special.jl

@@ -113,6 +113,8 @@ Random.seed!(1)
         Base.convert(::Type{TypeWithZero}, ::TypeWithoutZero) = TypeWithZero()
         Base.zero(::Type{<:Union{TypeWithoutZero, TypeWithZero}}) = TypeWithZero()
         LinearAlgebra.symmetric(::TypeWithoutZero, ::Symbol) = TypeWithoutZero()
+        LinearAlgebra.symmetric_type(::Type{TypeWithoutZero}) = TypeWithoutZero
+        Base.copy(A::TypeWithoutZero) = A
         Base.transpose(::TypeWithoutZero) = TypeWithoutZero()
         d  = fill(TypeWithoutZero(), 3)
         du = fill(TypeWithoutZero(), 2)

stdlib/LinearAlgebra/test/tridiag.jl

@@ -830,4 +830,12 @@ end
     @test axes(B) === (ax, ax)
 end
 
+@testset "Matrix conversion for non-numeric and undef" begin
+    T = Tridiagonal(fill(big(3), 3), Vector{BigInt}(undef, 4), fill(big(3), 3))
+    M = Matrix(T)
+    T[diagind(T)] .= 4
+    M[diagind(M)] .= 4
+    @test diag(T) == diag(M)
+end
+
 end # module TestTridiagonal