Higher-order functions in factorize to obtain structure

Author: jishnub

src/dense.jl

@@ -1541,38 +1541,8 @@ function factorize(A::AbstractMatrix{T}) where T
     m, n = size(A)
     if m == n
         if m == 1 return A[1] end
-        utri    = true
-        utri1   = true
-        herm    = true
-        sym     = true
-        for j = 1:n-1, i = j+1:m
-            if utri1
-                if A[i,j] != 0
-                    utri1 = i == j + 1
-                    utri = false
-                end
-            end
-            if sym
-                sym &= A[i,j] == A[j,i]
-            end
-            if herm
-                herm &= A[i,j] == conj(A[j,i])
-            end
-            if !(utri1|herm|sym) break end
-        end
-        ltri = true
-        ltri1 = true
-        for j = 3:n, i = 1:j-2
-            ltri1 &= A[i,j] == 0
-            if !ltri1 break end
-        end
+        utri, utri1, ltri, ltri1, sym, herm = getstructure(A)
         if ltri1
-            for i = 1:n-1
-                if A[i,i+1] != 0
-                    ltri &= false
-                    break
-                end
-            end
             if ltri
                 if utri
                     return Diagonal(A)
@@ -1610,6 +1580,54 @@ factorize(A::Adjoint)   =   adjoint(factorize(parent(A)))
 factorize(A::Transpose) = transpose(factorize(parent(A)))
 factorize(a::Number)    = a # same as how factorize behaves on Diagonal types
 
+function getstructure(A::StridedMatrix)
+    m, n = size(A)
+    if m == 1 return A[1] end
+    utri    = true
+    utri1   = true
+    herm    = true
+    sym     = true
+    for j = 1:n-1, i = j+1:m
+        if utri1
+            if A[i,j] != 0
+                utri1 = i == j + 1
+                utri = false
+            end
+        end
+        if sym
+            sym &= A[i,j] == A[j,i]
+        end
+        if herm
+            herm &= A[i,j] == conj(A[j,i])
+        end
+        if !(utri1|herm|sym) break end
+    end
+    ltri = true
+    ltri1 = true
+    for j = 3:n, i = 1:j-2
+        ltri1 &= A[i,j] == 0
+        if !ltri1 break end
+    end
+    if ltri1
+        for i = 1:n-1
+            if A[i,i+1] != 0
+                ltri &= false
+                break
+            end
+        end
+    end
+    return (utri, utri1, ltri, ltri1, sym, herm)
+end
+function getstructure(A::AbstractMatrix)
+    utri = istriu(A)
+    utri1 = istriu(A,-1)
+    ltri = istril(A)
+    ltri1 = istril(A,1)
+    sym = issymmetric(A)
+    herm = ishermitian(A)
+    return (utri, utri1, ltri, ltri1, sym, herm)
+end
+
 ## Moore-Penrose pseudoinverse
 
 """

test/dense.jl

@@ -4,6 +4,7 @@ module TestDense
 
 using Test, LinearAlgebra, Random
 using LinearAlgebra: BlasComplex, BlasFloat, BlasReal
+using Test: GenericArray
 
 const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")
 isdefined(Main, :FillArrays) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "FillArrays.jl"))
@@ -191,6 +192,8 @@ bimg  = randn(n,2)/2
         f = rand(eltya,n-2)
         A = diagm(0 => d)
         @test factorize(A) == Diagonal(d)
+        # test that the generic structure-evaluation method works
+        @test factorize(A) == factorize(GenericArray(A))
         A += diagm(-1 => e)
         @test factorize(A) == Bidiagonal(d,e,:L)
         A += diagm(-2 => f)