Higher-order functions in factorize to obtain structure (#1135)

Instead of looping over the entire array to check the matrix structure,
we may use higher-order querying functions like `istriu` and
`issymmetric`. The advantage of this is that matrices might have
optimized methods for these functions.

We still loop over the entire matrix for `StridedMatrix`es as before,
and obtain the structure in one-pass. The performance therefore isn't
impacted in this case.

An example with a sparse array wrapper:
```julia
julia> using LinearAlgebra

julia> struct MyMatrix{T,M<:AbstractMatrix{T}} <: AbstractMatrix{T}
           A::M
       end

julia> Base.size(M::MyMatrix) = size(M.A)

julia> Base.getindex(M::MyMatrix, i::Int, j::Int) = M.A[i, j]

julia> LinearAlgebra.istriu(M::MyMatrix, k::Integer=0) = istriu(M.A, k)

julia> LinearAlgebra.istril(M::MyMatrix, k::Integer=0) = istril(M.A, k)

julia> LinearAlgebra.issymmetric(M::MyMatrix) = issymmetric(M.A)

julia> LinearAlgebra.ishermitian(M::MyMatrix) = ishermitian(M.A)

julia> using SparseArrays

julia> S = sparse(1:4000, 1:4000, 1:4000);

julia> M = MyMatrix(S);

julia> @Btime factorize($M);
  178.231 ms (4 allocations: 31.34 KiB) # master
  22.165 ms (10 allocations: 94.04 KiB) # this PR
```

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,66 @@ 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
+_check_sym_herm(A) = (issymmetric(A), ishermitian(A))
+_check_sym_herm(A::AbstractMatrix{<:Real}) = (sym = issymmetric(A); (sym,sym))
+function getstructure(A::AbstractMatrix)
+    utri1 = istriu(A,-1)
+    # utri = istriu(A), but since we've already checked istriu(A,-1),
+    # we only need to check that the subdiagonal band is zero
+    utri = utri1 && iszero(diag(A,-1))
+    sym, herm = _check_sym_herm(A)
+    if sym || herm
+        # in either case, the lower and upper triangular halves have identical band structures
+        # in this case, istril(A,1) == istriu(A,-1) and istril(A) == istriu(A)
+        ltri1 = utri1
+        ltri = utri
+    else
+        ltri1 = istril(A,1)
+        # ltri = istril(A), but since we've already checked istril(A,1),
+        # we only need to check the superdiagonal band is zero
+        ltri = ltri1 && iszero(diag(A,1))
+    end
+    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)