add lq and tests

Author: Jutho

src/MatrixAlgebraKit.jl

@@ -31,6 +31,7 @@ export left_polar!, right_polar!
 export left_orth, right_orth, left_null, right_null
 export left_orth!, right_orth!, left_null!, right_null!
 
+export Native_HouseholderQR, Native_HouseholderLQ
 export LAPACK_HouseholderQR, LAPACK_HouseholderLQ, LAPACK_Simple, LAPACK_Expert,
     LAPACK_QRIteration, LAPACK_Bisection, LAPACK_MultipleRelativelyRobustRepresentations,
     LAPACK_DivideAndConquer, LAPACK_Jacobi

src/common/householder.jl

@@ -9,21 +9,21 @@ end
 Base.adjoint(H::Householder) = Householder(conj(H.β), H.v, H.r)
 
 function householder(x::AbstractVector, r::IndexRange = axes(x, 1), k = first(r))
-    i = findfirst(equalto(k), r)
+    i = findfirst(==(k), r)
     i == nothing && error("k = $k should be in the range r = $r")
     β, v, ν = _householder!(x[r], i)
     return Householder(β, v, r), ν
 end
 # Householder reflector h that zeros the elements A[r,col] (except for A[k,col]) upon lmul!(h,A)
 function householder(A::AbstractMatrix, r::IndexRange, col::Int, k = first(r))
-    i = findfirst(equalto(k), r)
+    i = findfirst(==(k), r)
     i == nothing && error("k = $k should be in the range r = $r")
     β, v, ν = _householder!(A[r, col], i)
     return Householder(β, v, r), ν
 end
 # Householder reflector that zeros the elements A[row,r] (except for A[row,k]) upon rmul!(A,h')
 function householder(A::AbstractMatrix, row::Int, r::IndexRange, k = first(r))
-    i = findfirst(equalto(k), r)
+    i = findfirst(==(k), r)
     i == nothing && error("k = $k should be in the range r = $r")
     β, v, ν = _householder!(conj!(A[row, r]), i)
     return Householder(β, v, r), ν

src/implementations/lq.jl

@@ -270,3 +270,85 @@ function _diagonal_lq!(
 end
 
 _diagonal_lq_null!(A::AbstractMatrix, N; positive::Bool = false) = N
+
+# Native logic
+# -------------
+function lq_full!(A::AbstractMatrix, LQ, alg::Native_HouseholderLQ)
+    check_input(lq_full!, A, LQ, alg)
+    L, Q = LQ
+    A === Q &&
+        throw(ArgumentError("inplace Q not supported with native LQ implementation"))
+    _native_lq!(A, L, Q; alg.kwargs...)
+    return L, Q
+end
+function lq_compact!(A::AbstractMatrix, LQ, alg::Native_HouseholderLQ)
+    check_input(lq_compact!, A, LQ, alg)
+    L, Q = LQ
+    A === Q &&
+        throw(ArgumentError("inplace Q not supported with native LQ implementation"))
+    _native_lq!(A, L, Q; alg.kwargs...)
+    return L, Q
+end
+function lq_null!(A::AbstractMatrix, N, alg::Native_HouseholderLQ)
+    check_input(lq_null!, A, N, alg)
+    _native_lq_null!(A, N; alg.kwargs...)
+    return N
+end
+
+function _native_lq!(
+        A::AbstractMatrix, L::AbstractMatrix, Q::AbstractMatrix;
+        positive::Bool = true # always true regardless of setting
+    )
+    m, n = size(A)
+    minmn = min(m, n)
+    @inbounds for i in 1:minmn
+        for j in 1:(i - 1)
+            L[i, j] = A[i, j]
+        end
+        β, v, L[i, i] = _householder!(conj!(view(A, i, i:n)), 1)
+        for j in (i + 1):size(L, 2)
+            L[i, j] = 0
+        end
+        H = Householder(conj(β), v, i:n)
+        rmul!(A, H; rows = (i + 1):m)
+        # A[i, i] == 1; store β instead
+        A[i, i] = β
+    end
+    # copy remaining rows for m > n
+    @inbounds for j in 1:size(L, 2)
+        for i in (minmn + 1):m
+            L[i, j] = A[i, j]
+        end
+    end
+    # build Q
+    one!(Q)
+    @inbounds for i in minmn:-1:1
+        β = A[i, i]
+        A[i, i] = 1
+        Hᴴ = Householder(β, view(A, i, i:n), i:n)
+        rmul!(Q, Hᴴ)
+    end
+    return L, Q
+end
+
+function _native_lq_null!(A::AbstractMatrix, Nᴴ::AbstractMatrix; positive::Bool = true)
+    m, n = size(A)
+    minmn = min(m, n)
+    @inbounds for i in 1:minmn
+        β, v, ν = _householder!(conj!(view(A, i, i:n)), 1)
+        H = Householder(conj(β), v, i:n)
+        rmul!(A, H; rows = (i + 1):m)
+        # A[i, i] == 1; store β instead
+        A[i, i] = β
+    end
+    # build Nᴴ
+    fill!(Nᴴ, zero(eltype(Nᴴ)))
+    one!(view(Nᴴ, 1:(n - minmn), (minmn + 1):n))
+    @inbounds for i in minmn:-1:1
+        β = A[i, i]
+        A[i, i] = 1
+        Hᴴ = Householder(β, view(A, i, i:n), i:n)
+        rmul!(Nᴴ, Hᴴ)
+    end
+    return Nᴴ
+end

src/implementations/qr.jl

@@ -350,7 +350,8 @@ function qr_null!(A::AbstractMatrix, N, alg::Native_HouseholderQR)
 end
 
 function _native_qr!(
-        A::AbstractMatrix, Q::AbstractMatrix, R::AbstractMatrix
+        A::AbstractMatrix, Q::AbstractMatrix, R::AbstractMatrix;
+        positive::Bool = true # always true regardless of setting
     )
     m, n = size(A)
     minmn = min(m, n)
@@ -367,14 +368,15 @@ function _native_qr!(
         # A[j,j] == 1; store β instead
         A[j, j] = β
     end
+    # copy remaining columns if m < n
     @inbounds for j in (minmn + 1):n
         for i in 1:size(R, 1)
             R[i, j] = A[i, j]
         end
     end
     # build Q
     one!(Q)
-    for j in minmn:-1:1
+    @inbounds for j in minmn:-1:1
         β = A[j, j]
         A[j, j] = 1
         Hᴴ = Householder(conj(β), view(A, j:m, j), j:m)
@@ -383,7 +385,7 @@ function _native_qr!(
     return Q, R
 end
 
-function _native_qr_null!(A::AbstractMatrix, N::AbstractMatrix)
+function _native_qr_null!(A::AbstractMatrix, N::AbstractMatrix; positive::Bool = true)
     m, n = size(A)
     minmn = min(m, n)
     @inbounds for j in 1:minmn
@@ -393,10 +395,10 @@ function _native_qr_null!(A::AbstractMatrix, N::AbstractMatrix)
         # A[j,j] == 1; store β instead
         A[j, j] = β
     end
-    # build Q
+    # build N
     fill!(N, zero(eltype(N)))
     one!(view(N, (minmn + 1):m, 1:(m - minmn)))
-    for j in minmn:-1:1
+    @inbounds for j in minmn:-1:1
         β = A[j, j]
         A[j, j] = 1
         Hᴴ = Householder(conj(β), view(A, j:m, j), j:m)

src/interface/decompositions.jl

@@ -10,14 +10,23 @@
 # QR, LQ, QL, RQ Decomposition
 # ----------------------------
 """
-    Native_HouseholderQR(; blocksize, positive = false, pivoted = false)
+    Native_HouseholderQR()
 
 Algorithm type to denote a native implementation for computing the QR decomposition of
-a matrix using Householder reflectors, .The diagonal elements of `R` will be non-negative
+a matrix using Householder reflectors. The diagonal elements of `R` will be non-negative
 by construction.
 """
 @algdef Native_HouseholderQR
 
+"""
+    Native_HouseholderLQ()
+
+Algorithm type to denote a native implementation for computing the LQ decomposition of
+a matrix using Householder reflectors. The diagonal elements of `L` will be non-negative
+by construction.
+"""
+@algdef Native_HouseholderLQ
+
 """
     LAPACK_HouseholderQR(; blocksize, positive = false, pivoted = false)
 

test/lq.jl

@@ -52,8 +52,8 @@ for T in (BLASFloats..., GenericFloats...), n in (37, m, 63)
             TestSuite.test_lq_algs(T, (m, n), LAPACK_LQ_ALGS)
         elseif T ∈ GenericFloats
             TestSuite.test_lq(T, (m, n); test_pivoted = false, test_blocksize = false)
-            GLA_LQ_ALGS = (LQViaTransposedQR(GLA_HouseholderQR()),)
-            TestSuite.test_lq_algs(T, (m, n), GLA_LQ_ALGS)
+            NATIVE_LQ_ALGS = (Native_HouseholderLQ(),)
+            TestSuite.test_lq_algs(T, (m, n), NATIVE_LQ_ALGS)
         end
         if m == n
             AT = Diagonal{T, Vector{T}}

test/qr.jl

@@ -51,8 +51,8 @@ for T in (BLASFloats..., GenericFloats...), n in (37, m, 63)
             TestSuite.test_qr_algs(T, (m, n), LAPACK_QR_ALGS)
         elseif T ∈ GenericFloats
             TestSuite.test_qr(T, (m, n); test_pivoted = false, test_blocksize = false)
-            GLA_QR_ALGS = (GLA_HouseholderQR(),)
-            TestSuite.test_qr_algs(T, (m, n), GLA_QR_ALGS)
+            NATIVE_QR_ALGS = (Native_HouseholderQR(),)
+            TestSuite.test_qr_algs(T, (m, n), NATIVE_QR_ALGS)
         end
         if m == n
             AT = Diagonal{T, Vector{T}}