rewrite diagm with Pairs (#24047)

* rewrite diagm with `Pair`s

Author: fredrikekre

NEWS.md

@@ -316,6 +316,8 @@ Library improvements
 
   * REPL Undo via Ctrl-/ and Ctrl-_
 
+  * `diagm` now accepts several diagonal index/vector `Pair`s ([#24047]).
+
   * New function `equalto(x)`, which returns a function that compares its argument to `x`
     using `isequal` ([#23812]).
 
@@ -467,10 +469,16 @@ Deprecated or removed
   * The tuple-of-types form of `cfunction`, `cfunction(f, returntype, (types...))`, has been deprecated
     in favor of the tuple-type form `cfunction(f, returntype, Tuple{types...})` ([#23066]).
 
+  * `diagm(v::AbstractVector, k::Integer=0)` has been deprecated in favor of
+    `diagm(k => v)` ([#24047]).
+
+  * `diagm(x::Number)` has been deprecated in favor of `fill(x, 1, 1)` ([#24047]).
+
   * `diagm(A::SparseMatrixCSC)` has been deprecated in favor of
     `spdiagm(sparsevec(A))` ([#23341]).
 
-  * `diagm(A::BitMatrix)` has been deprecated, use `diagm(vec(A))` instead ([#23373]).
+  * `diagm(A::BitMatrix)` has been deprecated, use `diagm(0 => vec(A))` or
+    `BitMatrix(Diagonal(vec(A)))` instead ([#23373], [#24047]).
 
   * `ℯ` (written as `\mscre<TAB>` or `\euler<TAB>`) is now the only (by default) exported
     name for Euler's number, and the type has changed from `Irrational{:e}` to

base/deprecated.jl

@@ -1641,7 +1641,7 @@ import .LinAlg: diagm
 @deprecate diagm(A::SparseMatrixCSC) sparse(Diagonal(sparsevec(A)))
 
 # PR #23373
-@deprecate diagm(A::BitMatrix) diagm(vec(A))
+@deprecate diagm(A::BitMatrix) BitMatrix(Diagonal(vec(A)))
 
 # PR 23341
 @eval GMP @deprecate gmp_version() version() false
@@ -1870,6 +1870,19 @@ end
     nothing
 end
 
+function diagm(v::BitVector)
+    depwarn(string("diagm(v::BitVector) is deprecated, use diagm(0 => v) or ",
+        "BitMatrix(Diagonal(v)) instead"), :diagm)
+    return BitMatrix(Diagonal(v))
+end
+function diagm(v::AbstractVector)
+    depwarn(string("diagm(v::AbstractVector) is deprecated, use diagm(0 => v) or ",
+        "Matrix(Diagonal(v)) instead"), :diagm)
+    return Matrix(Diagonal(v))
+end
+@deprecate diagm(v::AbstractVector, k::Integer) diagm(k => v)
+@deprecate diagm(x::Number) fill(x, 1, 1)
+
 # issue #20816
 @deprecate strwidth textwidth
 @deprecate charwidth textwidth

base/linalg/bitarray.jl

@@ -74,7 +74,7 @@ function tril(B::BitMatrix, k::Integer=0)
     A
 end
 
-## diag and related
+## diag
 
 function diag(B::BitMatrix)
     n = minimum(size(B))
@@ -85,15 +85,6 @@ function diag(B::BitMatrix)
     v
 end
 
-function diagm(v::BitVector)
-    n = length(v)
-    a = falses(n, n)
-    for i=1:n
-        a[i,i] = v[i]
-    end
-    a
-end
-
 ## norm and rank
 
 svd(A::BitMatrix) = svd(float(A))

base/linalg/dense.jl

@@ -246,7 +246,8 @@ diagind(A::AbstractMatrix, k::Integer=0) = diagind(size(A,1), size(A,2), k)
     diag(M, k::Integer=0)
 
 The `k`th diagonal of a matrix, as a vector.
-Use [`diagm`](@ref) to construct a diagonal matrix.
+
+See also: [`diagm`](@ref)
 
 # Examples
 ```jldoctest
@@ -265,29 +266,53 @@ julia> diag(A,1)
 diag(A::AbstractMatrix, k::Integer=0) = A[diagind(A,k)]
 
 """
-    diagm(v, k::Integer=0)
+    diagm(kv::Pair{<:Integer,<:AbstractVector}...)
+
+Construct a square matrix from `Pair`s of diagonals and vectors.
+Vector `kv.second` will be placed on the `kv.first` diagonal.
+`diagm` constructs a full matrix; if you want storage-efficient
+versions with fast arithmetic, see [`Diagonal`](@ref), [`Bidiagonal`](@ref)
+[`Tridiagonal`](@ref) and [`SymTridiagonal`](@ref).
 
-Construct a matrix by placing `v` on the `k`th diagonal. This constructs a full matrix; if
-you want a storage-efficient version with fast arithmetic, use [`Diagonal`](@ref) instead.
+See also: [`spdiagm`](@ref)
 
 # Examples
 ```jldoctest
-julia> diagm([1,2,3],1)
+julia> diagm(1 => [1,2,3])
 4×4 Array{Int64,2}:
  0  1  0  0
  0  0  2  0
  0  0  0  3
  0  0  0  0
+
+julia> diagm(1 => [1,2,3], -1 => [4,5])
+4×4 Array{Int64,2}:
+ 0  1  0  0
+ 4  0  2  0
+ 0  5  0  3
+ 0  0  0  0
 ```
 """
-function diagm(v::AbstractVector{T}, k::Integer=0) where T
-    n = length(v) + abs(k)
-    A = zeros(T,n,n)
-    A[diagind(A,k)] = v
-    A
+function diagm(kv::Pair{<:Integer,<:AbstractVector}...)
+    A = diagm_container(kv...)
+    for p in kv
+        inds = diagind(A, p.first)
+        for (i, val) in enumerate(p.second)
+            A[inds[i]] += val
+        end
+    end
+    return A
+end
+function diagm_container(kv::Pair{<:Integer,<:AbstractVector}...)
+    T = promote_type(map(x -> eltype(x.second), kv)...)
+    n = mapreduce(x -> length(x.second) + abs(x.first), max, kv)
+    return zeros(T, n, n)
+end
+function diagm_container(kv::Pair{<:Integer,<:BitVector}...)
+    n = mapreduce(x -> length(x.second) + abs(x.first), max, kv)
+    return falses(n, n)
 end
 
-diagm(x::Number) = (X = Matrix{typeof(x)}(1,1); X[1,1] = x; X)
 
 function trace(A::Matrix{T}) where T
     n = checksquare(A)

base/linalg/diagonal.jl

@@ -52,7 +52,7 @@ Diagonal{T}(V::AbstractVector) where {T} = Diagonal{T}(convert(AbstractVector{T}
 convert(::Type{Diagonal{T}}, D::Diagonal{T}) where {T} = D
 convert(::Type{Diagonal{T}}, D::Diagonal) where {T} = Diagonal{T}(convert(AbstractVector{T}, D.diag))
 convert(::Type{AbstractMatrix{T}}, D::Diagonal) where {T} = convert(Diagonal{T}, D)
-convert(::Type{Matrix}, D::Diagonal) = diagm(D.diag)
+convert(::Type{Matrix}, D::Diagonal) = diagm(0 => D.diag)
 convert(::Type{Array}, D::Diagonal) = convert(Matrix, D)
 full(D::Diagonal) = convert(Array, D)
 

base/linalg/generic.jl

@@ -739,13 +739,13 @@ of the [`eltype`](@ref) of `A`.
 julia> rank(eye(3))
 3
 
-julia> rank(diagm([1, 0, 2]))
+julia> rank(diagm(0 => [1, 0, 2]))
 2
 
-julia> rank(diagm([1, 0.001, 2]), 0.1)
+julia> rank(diagm(0 => [1, 0.001, 2]), 0.1)
 2
 
-julia> rank(diagm([1, 0.001, 2]), 0.00001)
+julia> rank(diagm(0 => [1, 0.001, 2]), 0.00001)
 3
 ```
 """

base/linalg/lapack.jl

@@ -1161,17 +1161,17 @@ of `A`. `fact` may be `E`, in which case `A` will be equilibrated and copied
 to `AF`; `F`, in which case `AF` and `ipiv` from a previous `LU` factorization
 are inputs; or `N`, in which case `A` will be copied to `AF` and then
 factored. If `fact = F`, `equed` may be `N`, meaning `A` has not been
-equilibrated; `R`, meaning `A` was multiplied by `diagm(R)` from the left;
-`C`, meaning `A` was multiplied by `diagm(C)` from the right; or `B`, meaning
-`A` was multiplied by `diagm(R)` from the left and `diagm(C)` from the right.
+equilibrated; `R`, meaning `A` was multiplied by `Diagonal(R)` from the left;
+`C`, meaning `A` was multiplied by `Diagonal(C)` from the right; or `B`, meaning
+`A` was multiplied by `Diagonal(R)` from the left and `Diagonal(C)` from the right.
 If `fact = F` and `equed = R` or `B` the elements of `R` must all be positive.
 If `fact = F` and `equed = C` or `B` the elements of `C` must all be positive.
 
 Returns the solution `X`; `equed`, which is an output if `fact` is not `N`,
 and describes the equilibration that was performed; `R`, the row equilibration
 diagonal; `C`, the column equilibration diagonal; `B`, which may be overwritten
-with its equilibrated form `diagm(R)*B` (if `trans = N` and `equed = R,B`) or
-`diagm(C)*B` (if `trans = T,C` and `equed = C,B`); `rcond`, the reciprocal
+with its equilibrated form `Diagonal(R)*B` (if `trans = N` and `equed = R,B`) or
+`Diagonal(C)*B` (if `trans = T,C` and `equed = C,B`); `rcond`, the reciprocal
 condition number of `A` after equilbrating; `ferr`, the forward error bound for
 each solution vector in `X`; `berr`, the forward error bound for each solution
 vector in `X`; and `work`, the reciprocal pivot growth factor.

base/linalg/svd.jl

@@ -32,7 +32,7 @@ end
 Compute the singular value decomposition (SVD) of `A` and return an `SVD` object.
 
 `U`, `S`, `V` and `Vt` can be obtained from the factorization `F` with `F[:U]`,
-`F[:S]`, `F[:V]` and `F[:Vt]`, such that `A = U*diagm(S)*Vt`.
+`F[:S]`, `F[:V]` and `F[:Vt]`, such that `A = U*Diagonal(S)*Vt`.
 The algorithm produces `Vt` and hence `Vt` is more efficient to extract than `V`.
 The singular values in `S` are sorted in descending order.
 
@@ -52,7 +52,7 @@ julia> A = [1. 0. 0. 0. 2.; 0. 0. 3. 0. 0.; 0. 0. 0. 0. 0.; 0. 2. 0. 0. 0.]
 julia> F = svdfact(A)
 Base.LinAlg.SVD{Float64,Float64,Array{Float64,2}}([0.0 1.0 0.0 0.0; 1.0 0.0 0.0 0.0; 0.0 0.0 0.0 -1.0; 0.0 0.0 1.0 0.0], [3.0, 2.23607, 2.0, 0.0], [-0.0 0.0 … -0.0 0.0; 0.447214 0.0 … 0.0 0.894427; -0.0 1.0 … -0.0 0.0; 0.0 0.0 … 1.0 0.0])
 
-julia> F[:U] * diagm(F[:S]) * F[:Vt]
+julia> F[:U] * Diagonal(F[:S]) * F[:Vt]
 4×5 Array{Float64,2}:
  1.0  0.0  0.0  0.0  2.0
  0.0  0.0  3.0  0.0  0.0
@@ -71,7 +71,7 @@ svdfact(x::Integer; thin::Bool=true) = svdfact(float(x), thin=thin)
     svd(A; thin::Bool=true) -> U, S, V
 
 Computes the SVD of `A`, returning `U`, vector `S`, and `V` such that
-`A == U*diagm(S)*V'`. The singular values in `S` are sorted in descending order.
+`A == U*Diagonal(S)*V'`. The singular values in `S` are sorted in descending order.
 
 If `thin=true` (default), a thin SVD is returned. For a ``M \\times N`` matrix
 `A`, `U` is ``M \\times M`` for a full SVD (`thin=false`) and
@@ -93,7 +93,7 @@ julia> A = [1. 0. 0. 0. 2.; 0. 0. 3. 0. 0.; 0. 0. 0. 0. 0.; 0. 2. 0. 0. 0.]
 julia> U, S, V = svd(A)
 ([0.0 1.0 0.0 0.0; 1.0 0.0 0.0 0.0; 0.0 0.0 0.0 -1.0; 0.0 0.0 1.0 0.0], [3.0, 2.23607, 2.0, 0.0], [-0.0 0.447214 -0.0 0.0; 0.0 0.0 1.0 0.0; … ; -0.0 0.0 -0.0 1.0; 0.0 0.894427 0.0 0.0])
 
-julia> U*diagm(S)*V'
+julia> U*Diagonal(S)*V'
 4×5 Array{Float64,2}:
  1.0  0.0  0.0  0.0  2.0
  0.0  0.0  3.0  0.0  0.0
@@ -266,17 +266,17 @@ function getindex(obj::GeneralizedSVD{T}, d::Symbol) where T
     elseif d == :D1
         m = size(obj.U, 1)
         if m - obj.k - obj.l >= 0
-            return [eye(T, obj.k) zeros(T, obj.k, obj.l); zeros(T, obj.l, obj.k) diagm(obj.a[obj.k + 1:obj.k + obj.l]); zeros(T, m - obj.k - obj.l, obj.k + obj.l)]
+            return [eye(T, obj.k) zeros(T, obj.k, obj.l); zeros(T, obj.l, obj.k) Diagonal(obj.a[obj.k + 1:obj.k + obj.l]); zeros(T, m - obj.k - obj.l, obj.k + obj.l)]
         else
-            return [eye(T, m, obj.k) [zeros(T, obj.k, m - obj.k); diagm(obj.a[obj.k + 1:m])] zeros(T, m, obj.k + obj.l - m)]
+            return [eye(T, m, obj.k) [zeros(T, obj.k, m - obj.k); Diagonal(obj.a[obj.k + 1:m])] zeros(T, m, obj.k + obj.l - m)]
         end
     elseif d == :D2
         m = size(obj.U, 1)
         p = size(obj.V, 1)
         if m - obj.k - obj.l >= 0
-            return [zeros(T, obj.l, obj.k) diagm(obj.b[obj.k + 1:obj.k + obj.l]); zeros(T, p - obj.l, obj.k + obj.l)]
+            return [zeros(T, obj.l, obj.k) Diagonal(obj.b[obj.k + 1:obj.k + obj.l]); zeros(T, p - obj.l, obj.k + obj.l)]
         else
-            return [zeros(T, p, obj.k) [diagm(obj.b[obj.k + 1:m]); zeros(T, obj.k + p - m, m - obj.k)] [zeros(T, m - obj.k, obj.k + obj.l - m); eye(T, obj.k + p - m, obj.k + obj.l - m)]]
+            return [zeros(T, p, obj.k) [Diagonal(obj.b[obj.k + 1:m]); zeros(T, obj.k + p - m, m - obj.k)] [zeros(T, m - obj.k, obj.k + obj.l - m); eye(T, obj.k + p - m, obj.k + obj.l - m)]]
         end
     elseif d == :R
         return obj.R

base/linalg/uniformscaling.jl

@@ -220,7 +220,7 @@ function isapprox(J::UniformScaling,A::AbstractMatrix;
                   rtol::Real=rtoldefault(promote_leaf_eltypes(A),eltype(J),atol),
                   nans::Bool=false, norm::Function=vecnorm)
     n = checksquare(A)
-    Jnorm = norm === vecnorm ? abs(J.λ)*sqrt(n) : (norm === Base.norm ? abs(J.λ) : norm(diagm(fill(J.λ, n))))
+    Jnorm = norm === vecnorm ? abs(J.λ)*sqrt(n) : (norm === Base.norm ? abs(J.λ) : norm(Diagonal(fill(J.λ, n))))
     return norm(A - J) <= max(atol, rtol*max(norm(A), Jnorm))
 end
 isapprox(A::AbstractMatrix,J::UniformScaling;kwargs...) = isapprox(J,A;kwargs...)

test/arrayops.jl

@@ -232,7 +232,7 @@ end
     tmp[rng...] = A[rng...]
     @test  tmp == cat(3,zeros(Int,2,3),[0 0 0; 0 47 52],zeros(Int,2,3),[0 0 0; 0 127 132])
 
-    @test cat([1,2],1,2,3.,4.,5.) == diagm([1,2,3.,4.,5.])
+    @test cat([1,2],1,2,3.,4.,5.) == diagm(0 => [1,2,3.,4.,5.])
     blk = [1 2;3 4]
     tmp = cat([1,3],blk,blk)
     @test tmp[1:2,1:2,1] == blk
@@ -2044,7 +2044,7 @@ end # module AutoRetType
 @testset "concatenations of dense matrices/vectors yield dense matrices/vectors" begin
     N = 4
     densevec = ones(N)
-    densemat = diagm(ones(N))
+    densemat = diagm(0 => ones(N))
     # Test that concatenations of homogeneous pairs of either dense matrices or dense vectors
     # (i.e., Matrix-Matrix concatenations, and Vector-Vector concatenations) yield dense arrays
     for densearray in (densevec, densemat)