implement fillband! and call it from fillstored!(::AbstractTriangular)

Author: fredrikekre

base/linalg/dense.jl

@@ -242,6 +242,26 @@ function tril!(M::AbstractMatrix, k::Integer)
 end
 tril(M::Matrix, k::Integer) = tril!(copy(M), k)
 
+"""
+    fillband!(A::AbstractMatrix, x, l, u)
+
+Fill the band between diagonals `l` and `u` with the value `x`.
+"""
+function fillband!(A::AbstractMatrix{T}, x, l, u) where T
+    m, n = size(A)
+    xT = convert(T, x)
+    idx = 1
+    for j in 0:n-1
+        iu = min(max(0, j-u), m) + idx
+        il = min(max(0, j-l), m-1) + idx
+        for i in iu:il
+            A[i] = xT
+        end
+        idx += m
+    end
+    return A
+end
+
 function diagind(m::Integer, n::Integer, k::Integer=0)
     if !(-m <= k <= n)
         throw(ArgumentError(string("requested diagonal, $k, must be at least $(-m) and ",

base/linalg/triangular.jl

@@ -418,7 +418,10 @@ end
 scale!(A::Union{UpperTriangular,LowerTriangular}, c::Number) = scale!(A,A,c)
 scale!(c::Number, A::Union{UpperTriangular,LowerTriangular}) = scale!(A,c)
 
-fillstored!(A::AbstractTriangular, x) = (fill!(A.data, x); A)
+fillstored!(A::LowerTriangular, x)     = (fillband!(A.data, x, size(A,1), 0); A)
+fillstored!(A::UnitLowerTriangular, x) = (fillband!(A.data, x, size(A,1), -1); A)
+fillstored!(A::UpperTriangular, x)     = (fillband!(A.data, x, 0, size(A,2)); A)
+fillstored!(A::UnitUpperTriangular, x) = (fillband!(A.data, x, 1, size(A,2)); A)
 
 # Binary operations
 +(A::UpperTriangular, B::UpperTriangular) = UpperTriangular(A.data + B.data)