Make fillstored! public (#1333)

`fillstored!` offers a way to fill only the stored indices of a
structured matrix without having to populate the parent. E.g.:
```julia
julia> U = UpperTriangular(zeros(4,4))
4×4 UpperTriangular{Float64, Matrix{Float64}}:
 0.0  0.0  0.0  0.0
  ⋅   0.0  0.0  0.0
  ⋅    ⋅   0.0  0.0
  ⋅    ⋅    ⋅   0.0

julia> LinearAlgebra.fillstored!(U, 2)
4×4 UpperTriangular{Float64, Matrix{Float64}}:
 2.0  2.0  2.0  2.0
  ⋅   2.0  2.0  2.0
  ⋅    ⋅   2.0  2.0
  ⋅    ⋅    ⋅   2.0
```
This seems like a useful function that should be public. This came up on
discourse:
https://discourse.julialang.org/t/how-to-set-all-elements-in-a-lower-triangular-matrix/128603

Author: jishnub

src/LinearAlgebra.jl

@@ -179,7 +179,8 @@ public AbstractTriangular,
         symmetric,
         symmetric_type,
         zeroslike,
-        matprod_dest
+        matprod_dest,
+        fillstored!
 
 const BlasFloat = Union{Float64,Float32,ComplexF64,ComplexF32}
 const BlasReal = Union{Float64,Float32}

src/special.jl

@@ -320,10 +320,33 @@ _diag_or_value(A::Diagonal) = A.diag
 _diag_or_value(A::UniformScaling) = A.λ
 
 # fill[stored]! methods
+"""
+    fillstored!(A::AbstractMatrix, x)
+
+Fill only the stored indices of a structured matrix `A` with the value `x`.
+
+# Example
+```jldoctest
+julia> A = Tridiagonal(zeros(2), zeros(3), zeros(2))
+3×3 Tridiagonal{Float64, Vector{Float64}}:
+ 0.0  0.0   ⋅
+ 0.0  0.0  0.0
+  ⋅   0.0  0.0
+
+julia> LinearAlgebra.fillstored!(A, 2)
+3×3 Tridiagonal{Float64, Vector{Float64}}:
+ 2.0  2.0   ⋅
+ 2.0  2.0  2.0
+  ⋅   2.0  2.0
+```
+"""
 fillstored!(A::Diagonal, x) = (fill!(A.diag, x); A)
 fillstored!(A::Bidiagonal, x) = (fill!(A.dv, x); fill!(A.ev, x); A)
 fillstored!(A::Tridiagonal, x) = (fill!(A.dl, x); fill!(A.d, x); fill!(A.du, x); A)
-fillstored!(A::SymTridiagonal, x) = (fill!(A.dv, x); fill!(A.ev, x); A)
+function fillstored!(A::SymTridiagonal, x)
+    issymmetric(x) || throw(ArgumentError("cannot set a diagonal entry of a SymTridiagonal to an asymmetric value"))
+    (fill!(A.dv, x); fill!(A.ev, x); A)
+end
 
 _small_enough(A::Union{Diagonal, Bidiagonal}) = size(A, 1) <= 1
 _small_enough(A::Tridiagonal) = size(A, 1) <= 2

src/symmetric.jl

@@ -401,13 +401,11 @@ fill!(A::HermOrSym, x) = fillstored!(A, x)
 function fillstored!(A::HermOrSym{T}, x) where T
     xT = convert(T, x)
     if isa(A, Hermitian)
-        isreal(xT) || throw(ArgumentError("cannot fill Hermitian matrix with a nonreal value"))
-    end
-    if A.uplo == 'U'
-        fillband!(A.data, xT, 0, size(A,2)-1)
-    else # A.uplo == 'L'
-        fillband!(A.data, xT, 1-size(A,1), 0)
+        ishermitian(xT) || throw(ArgumentError("cannot fill Hermitian matrix with a non-hermitian value"))
+    elseif isa(A, Symmetric)
+        issymmetric(xT) || throw(ArgumentError("cannot fill Symmetric matrix with an asymmetric value"))
     end
+    applytri(A -> fillstored!(A, xT), A)
     return A
 end
 

test/special.jl

@@ -835,4 +835,28 @@ end
     end
 end
 
+@testset "fillstored!" begin
+    dv, ev = zeros(4), zeros(3)
+    D = Diagonal(dv)
+    LinearAlgebra.fillstored!(D, 2)
+    @test D == diagm(fill(2, length(dv)))
+
+    dv .= 0
+    B = Bidiagonal(dv, ev, :U)
+    LinearAlgebra.fillstored!(B, 2)
+    @test B == diagm(0=>fill(2, length(dv)), 1=>fill(2, length(ev)))
+
+    dv .= 0
+    ev .= 0
+    T = Tridiagonal(ev, dv, ev)
+    LinearAlgebra.fillstored!(T, 2)
+    @test T == diagm(-1=>fill(2, length(ev)), 0=>fill(2, length(dv)), 1=>fill(2, length(ev)))
+
+    dv .= 0
+    ev .= 0
+    ST = SymTridiagonal(dv, ev)
+    LinearAlgebra.fillstored!(ST, 2)
+    @test ST == diagm(-1=>fill(2, length(ev)), 0=>fill(2, length(dv)), 1=>fill(2, length(ev)))
+end
+
 end # module TestSpecial

test/symmetric.jl

@@ -1165,4 +1165,13 @@ end
     end
 end
 
+@testset "fillstored!" begin
+    A = zeros(4,4)
+    for T in (Symmetric, Hermitian)
+        A .= 0
+        LinearAlgebra.fillstored!(T(A), 2)
+        @test all(==(2), T(A))
+    end
+end
+
 end # module TestSymmetric