From 2f9a91576f6e456329a6284f05244160681abb34 Mon Sep 17 00:00:00 2001
From: Jutho Haegeman <jutho.haegeman@ugent.be>
Date: Sat, 1 Nov 2025 00:50:32 +0100
Subject: [PATCH 1/3] native_qr

---
 src/MatrixAlgebraKit.jl         |   1 +
 src/common/householder.jl       | 142 ++++++++++++++++++++++++++++++++
 src/implementations/qr.jl       |  87 ++++++++++++++++++-
 src/interface/decompositions.jl |   9 ++
 4 files changed, 235 insertions(+), 4 deletions(-)
 create mode 100644 src/common/householder.jl

diff --git a/src/MatrixAlgebraKit.jl b/src/MatrixAlgebraKit.jl
index 5211476f..6d6a8065 100644
--- a/src/MatrixAlgebraKit.jl
+++ b/src/MatrixAlgebraKit.jl
@@ -69,6 +69,7 @@ export notrunc, truncrank, trunctol, truncerror, truncfilter
 end
 
 include("common/defaults.jl")
+include("common/householder.jl")
 include("common/initialization.jl")
 include("common/pullbacks.jl")
 include("common/safemethods.jl")
diff --git a/src/common/householder.jl b/src/common/householder.jl
new file mode 100644
index 00000000..6b1a1217
--- /dev/null
+++ b/src/common/householder.jl
@@ -0,0 +1,142 @@
+const IndexRange{T <: Integer} = Base.AbstractRange{T}
+
+# Elementary Householder reflection
+struct Householder{T, V <: AbstractVector, R <: IndexRange}
+    β::T
+    v::V
+    r::R
+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 == 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 == 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 == nothing && error("k = $k should be in the range r = $r")
+    β, v, ν = _householder!(conj!(A[row, r]), i)
+    return Householder(β, v, r), ν
+end
+
+# generate Householder vector based on vector v, such that applying the reflection
+# to v yields a vector with single non-zero element on position i, whose value is
+# positive and thus equal to norm(v)
+function _householder!(v::AbstractVector{T}, i::Int = 1) where {T}
+    β::T = zero(T)
+    @inbounds begin
+        σ = abs2(zero(T))
+        @simd for k in 1:(i - 1)
+            σ += abs2(v[k])
+        end
+        @simd for k in (i + 1):length(v)
+            σ += abs2(v[k])
+        end
+        vi = v[i]
+        ν = sqrt(abs2(vi) + σ)
+
+        if σ == 0 && vi == ν
+            β = zero(vi)
+        else
+            if real(vi) < 0
+                vi = vi - ν
+            else
+                vi = ((vi - conj(vi)) * ν - σ) / (conj(vi) + ν)
+            end
+            @simd for k in 1:(i - 1)
+                v[k] /= vi
+            end
+            v[i] = 1
+            @simd for k in (i + 1):length(v)
+                v[k] /= vi
+            end
+            β = -conj(vi) / (ν)
+        end
+    end
+    return β, v, ν
+end
+
+function LinearAlgebra.lmul!(H::Householder, x::AbstractVector)
+    v = H.v
+    r = H.r
+    β = H.β
+    β == 0 && return x
+    @inbounds begin
+        μ = conj(zero(v[1])) * zero(x[r[1]])
+        i = 1
+        @simd for j in r
+            μ += conj(v[i]) * x[j]
+            i += 1
+        end
+        μ *= β
+        i = 1
+        @simd for j in H.r
+            x[j] -= μ * v[i]
+            i += 1
+        end
+    end
+    return x
+end
+function LinearAlgebra.lmul!(H::Householder, A::AbstractMatrix; cols = axes(A, 2))
+    v = H.v
+    r = H.r
+    β = H.β
+    β == 0 && return A
+    @inbounds begin
+        for k in cols
+            μ = conj(zero(v[1])) * zero(A[r[1], k])
+            i = 1
+            @simd for j in r
+                μ += conj(v[i]) * A[j, k]
+                i += 1
+            end
+            μ *= β
+            i = 1
+            @simd for j in H.r
+                A[j, k] -= μ * v[i]
+                i += 1
+            end
+        end
+    end
+    return A
+end
+function LinearAlgebra.rmul!(A::AbstractMatrix, H::Householder; rows = axes(A, 1))
+    v = H.v
+    r = H.r
+    β = H.β
+    β == 0 && return A
+    w = similar(A, length(rows))
+    fill!(w, 0)
+    all(in(axes(A, 2)), r) || error("Householder range r = $r not compatible with matrix A of size $(size(A))")
+    @inbounds begin
+        l = 1
+        for k in r
+            j = 1
+            @simd for i in rows
+                w[j] += A[i, k] * v[l]
+                j += 1
+            end
+            l += 1
+        end
+        l = 1
+        for k in r
+            j = 1
+            @simd for i in rows
+                A[i, k] -= β * w[j] * conj(v[l])
+                j += 1
+            end
+            l += 1
+        end
+    end
+    return A
+end
diff --git a/src/implementations/qr.jl b/src/implementations/qr.jl
index a3c9d1f5..9e87779b 100644
--- a/src/implementations/qr.jl
+++ b/src/implementations/qr.jl
@@ -233,10 +233,9 @@ end
 
 _diagonal_qr_null!(A::AbstractMatrix, N; positive::Bool = false) = N
 
-### GPU logic
-# placed here to avoid code duplication since much of the logic is replicable across
-# CUDA and AMDGPU
-###
+# GPU logic
+# --------------
+# placed here to avoid code duplication since much of the logic is replicable across CUDA and AMDGPU
 function MatrixAlgebraKit.qr_full!(
         A::AbstractMatrix, QR, alg::Union{CUSOLVER_HouseholderQR, ROCSOLVER_HouseholderQR}
     )
@@ -321,3 +320,83 @@ function _gpu_qr_null!(
     N = _gpu_unmqr!('L', 'N', A, τ, N)
     return N
 end
+
+# Native logic
+# --------------
+function qr_full!(A::AbstractMatrix, QR, alg::Native_HouseholderQR)
+    check_input(qr_full!, A, QR, alg)
+    Q, R = QR
+    A === Q &&
+        throw(ArgumentError("inplace Q not supported with native QR implementation"))
+    _native_qr!(A, Q, R; alg.kwargs...)
+    return Q, R
+end
+function qr_compact!(A::AbstractMatrix, QR, alg::Native_HouseholderQR)
+    check_input(qr_compact!, A, QR, alg)
+    Q, R = QR
+    A === Q &&
+        throw(ArgumentError("inplace Q not supported with native QR implementation"))
+    _native_qr!(A, Q, R; alg.kwargs...)
+    return Q, R
+end
+function qr_null!(A::AbstractMatrix, N, alg::Native_HouseholderQR)
+    check_input(qr_null!, A, N, alg)
+    _native_qr_null!(A, N; alg.kwargs...)
+    return N
+end
+
+function _native_qr!(
+        A::AbstractMatrix, Q::AbstractMatrix, R::AbstractMatrix
+    )
+    m, n = size(A)
+    minmn = min(m, n)
+    @inbounds for j in 1:minmn
+        for i in 1:(j - 1)
+            R[i, j] = A[i, j]
+        end
+        β, v, R[j, j] = _householder!(view(A, j:m, j), 1)
+        for i in (j + 1):size(R, 1)
+            R[i, j] = 0
+        end
+        H = Householder(β, v, j:m)
+        lmul!(H, A; cols = (j + 1):n)
+        # A[j,j] == 1; store β instead
+        A[j, j] = β
+    end
+    @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
+        β = A[j, j]
+        A[j, j] = 1
+        Hᴴ = Householder(conj(β), view(A, j:m, j), j:m)
+        lmul!(Hᴴ, Q)
+    end
+    return Q, R
+end
+
+function _native_qr_null!(A::AbstractMatrix, N::AbstractMatrix)
+    m, n = size(A)
+    minmn = min(m, n)
+    @inbounds for j in 1:minmn
+        β, v, ν = _householder!(view(A, j:m, j), 1)
+        H = Householder(β, v, j:m)
+        lmul!(H, A; cols = (j + 1):n)
+        # A[j,j] == 1; store β instead
+        A[j, j] = β
+    end
+    # build Q
+    fill!(N, zero(eltype(N)))
+    one!(view(N, (minmn + 1):m, 1:(m - minmn)))
+    for j in minmn:-1:1
+        β = A[j, j]
+        A[j, j] = 1
+        Hᴴ = Householder(conj(β), view(A, j:m, j), j:m)
+        lmul!(Hᴴ, N)
+    end
+    return N
+end
diff --git a/src/interface/decompositions.jl b/src/interface/decompositions.jl
index bdda6612..f311caab 100644
--- a/src/interface/decompositions.jl
+++ b/src/interface/decompositions.jl
@@ -9,6 +9,15 @@
 
 # QR, LQ, QL, RQ Decomposition
 # ----------------------------
+"""
+    Native_HouseholderQR(; blocksize, positive = false, pivoted = false)
+
+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
+by construction.
+"""
+@algdef Native_HouseholderQR
+
 """
     LAPACK_HouseholderQR(; blocksize, positive = false, pivoted = false)
 

From 75e0d03f15bdba920efda2ef1fd07e3ddc36f4f1 Mon Sep 17 00:00:00 2001
From: Jutho Haegeman <jutho.haegeman@ugent.be>
Date: Sun, 2 Nov 2025 09:31:32 +0100
Subject: [PATCH 2/3] add lq and tests

---
 src/MatrixAlgebraKit.jl         |  1 +
 src/common/householder.jl       |  6 +--
 src/implementations/lq.jl       | 82 +++++++++++++++++++++++++++++++++
 src/implementations/qr.jl       | 12 +++--
 src/interface/decompositions.jl | 13 +++++-
 test/lq.jl                      | 57 +++++++++++++++++++++--
 test/qr.jl                      | 59 ++++++++++++++++++++++--
 7 files changed, 212 insertions(+), 18 deletions(-)

diff --git a/src/MatrixAlgebraKit.jl b/src/MatrixAlgebraKit.jl
index 6d6a8065..852cc3d8 100644
--- a/src/MatrixAlgebraKit.jl
+++ b/src/MatrixAlgebraKit.jl
@@ -30,6 +30,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
diff --git a/src/common/householder.jl b/src/common/householder.jl
index 6b1a1217..0aac21bc 100644
--- a/src/common/householder.jl
+++ b/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), ν
diff --git a/src/implementations/lq.jl b/src/implementations/lq.jl
index a1648d3e..33d4f132 100644
--- a/src/implementations/lq.jl
+++ b/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
diff --git a/src/implementations/qr.jl b/src/implementations/qr.jl
index 9e87779b..a8cd298a 100644
--- a/src/implementations/qr.jl
+++ b/src/implementations/qr.jl
@@ -346,7 +346,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)
@@ -363,6 +364,7 @@ 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]
@@ -370,7 +372,7 @@ function _native_qr!(
     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)
@@ -379,7 +381,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
@@ -389,10 +391,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)
diff --git a/src/interface/decompositions.jl b/src/interface/decompositions.jl
index f311caab..7c67118a 100644
--- a/src/interface/decompositions.jl
+++ b/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)
 
diff --git a/test/lq.jl b/test/lq.jl
index 2c8dfefe..31ce2696 100644
--- a/test/lq.jl
+++ b/test/lq.jl
@@ -5,9 +5,10 @@ using StableRNGs
 using LinearAlgebra: diag, I, Diagonal
 using MatrixAlgebraKit: LQViaTransposedQR, LAPACK_HouseholderQR
 
-eltypes = (Float32, Float64, ComplexF32, ComplexF64)
+lapack_eltypes = (Float32, Float64, ComplexF32, ComplexF64)
+native_eltypes = (lapack_eltypes..., BigFloat, Complex{BigFloat})
 
-@testset "lq_compact! for T = $T" for T in eltypes
+@testset "lq_compact! for T = $T" for T in lapack_eltypes
     rng = StableRNG(123)
     m = 54
     for n in (37, m, 63)
@@ -114,7 +115,7 @@ eltypes = (Float32, Float64, ComplexF32, ComplexF64)
     end
 end
 
-@testset "lq_full! for T = $T" for T in eltypes
+@testset "lq_full! for T = $T" for T in lapack_eltypes
     rng = StableRNG(123)
     m = 54
     for n in (37, m, 63)
@@ -208,7 +209,7 @@ end
     end
 end
 
-@testset "lq_compact, lq_full and lq_null for Diagonal{$T}" for T in eltypes
+@testset "lq_compact, lq_full and lq_null for Diagonal{$T}" for T in native_eltypes
     rng = StableRNG(123)
     atol = eps(real(T))^(3 / 4)
     for m in (54, 0)
@@ -252,3 +253,51 @@ end
         @test N isa AbstractMatrix{T} && size(N) == (0, m)
     end
 end
+
+@testset "native lq_compact! for T = $T" for T in native_eltypes
+    rng = StableRNG(123)
+    m = 54
+    for n in (37, m, 63)
+        minmn = min(m, n)
+        A = randn(rng, T, m, n)
+        L, Q = @constinferred lq_compact(A, Native_HouseholderLQ())
+        @test L isa Matrix{T} && size(L) == (m, minmn)
+        @test Q isa Matrix{T} && size(Q) == (minmn, n)
+        @test L * Q ≈ A
+        @test all(>=(zero(real(T))), real(diag(L)))
+        @test isisometric(Q; side = :right)
+        Nᴴ = @constinferred lq_null(A, Native_HouseholderLQ())
+        @test Nᴴ isa Matrix{T} && size(Nᴴ) == (n - minmn, n)
+        @test maximum(abs, A * Nᴴ') < eps(real(T))^(2 / 3)
+        @test isisometric(Nᴴ; side = :right)
+
+        Ac = similar(A)
+        L2, Q2 = @constinferred lq_compact!(copy!(Ac, A), (L, Q), Native_HouseholderLQ())
+        @test L2 === L
+        @test Q2 === Q
+        Nᴴ2 = @constinferred lq_null!(copy!(Ac, A), Nᴴ, Native_HouseholderLQ())
+        @test Nᴴ2 === Nᴴ
+    end
+end
+
+@testset "lq_full! for T = $T" for T in native_eltypes
+    rng = StableRNG(123)
+    m = 54
+    for n in (37, m, 63)
+        minmn = min(m, n)
+        A = randn(rng, T, m, n)
+        L, Q = lq_full(A, Native_HouseholderLQ())
+        @test L isa Matrix{T} && size(L) == (m, n)
+        @test Q isa Matrix{T} && size(Q) == (n, n)
+        @test L * Q ≈ A
+        @test all(>=(zero(real(T))), real(diag(L)))
+        @test isunitary(Q)
+
+        Ac = similar(A)
+        L2, Q2 = @constinferred lq_full!(copy!(Ac, A), (L, Q), Native_HouseholderLQ())
+        @test L2 === L
+        @test Q2 === Q
+        @test L * Q ≈ A
+        @test isunitary(Q)
+    end
+end
diff --git a/test/qr.jl b/test/qr.jl
index 826c320b..6e88dfc3 100644
--- a/test/qr.jl
+++ b/test/qr.jl
@@ -4,9 +4,10 @@ using TestExtras
 using StableRNGs
 using LinearAlgebra: diag, I, Diagonal
 
-eltypes = (Float32, Float64, ComplexF32, ComplexF64)
+lapack_eltypes = (Float32, Float64, ComplexF32, ComplexF64)
+native_eltypes = (lapack_eltypes..., BigFloat, Complex{BigFloat})
 
-@testset "qr_compact! and qr_null! for T = $T" for T in eltypes
+@testset "qr_compact! and qr_null! for T = $T" for T in lapack_eltypes
     rng = StableRNG(123)
     m = 54
     for n in (37, m, 63)
@@ -99,7 +100,7 @@ eltypes = (Float32, Float64, ComplexF32, ComplexF64)
     end
 end
 
-@testset "qr_full! for T = $T" for T in eltypes
+@testset "qr_full! for T = $T" for T in lapack_eltypes
     rng = StableRNG(123)
     m = 54
     for n in (37, m, 63)
@@ -176,7 +177,7 @@ end
     end
 end
 
-@testset "qr_compact, qr_full and qr_null for Diagonal{$T}" for T in eltypes
+@testset "qr_compact, qr_full and qr_null for Diagonal{$T}" for T in native_eltypes
     rng = StableRNG(123)
     atol = eps(real(T))^(3 / 4)
     for m in (54, 0)
@@ -220,3 +221,53 @@ end
         @test N isa AbstractMatrix{T} && size(N) == (m, 0)
     end
 end
+
+
+@testset "native qr_compact! and qr_null! for T = $T" for T in native_eltypes
+    rng = StableRNG(123)
+    m = 54
+    for n in (37, m, 63)
+        minmn = min(m, n)
+        A = randn(rng, T, m, n)
+        Q, R = @constinferred qr_compact(A, Native_HouseholderQR())
+        @test Q isa Matrix{T} && size(Q) == (m, minmn)
+        @test R isa Matrix{T} && size(R) == (minmn, n)
+        @test Q * R ≈ A
+        @test all(>=(zero(real(T))), real(diag(R)))
+        N = @constinferred qr_null(A, Native_HouseholderQR())
+        @test N isa Matrix{T} && size(N) == (m, m - minmn)
+        @test isisometric(Q)
+        @test maximum(abs, A' * N) < eps(real(T))^(2 / 3)
+        @test isisometric(N)
+
+        Ac = similar(A)
+        Q2, R2 = @constinferred qr_compact!(copy!(Ac, A), (Q, R), Native_HouseholderQR())
+        @test Q2 === Q
+        @test R2 === R
+        N2 = @constinferred qr_null!(copy!(Ac, A), N, Native_HouseholderQR())
+        @test N2 === N
+    end
+end
+
+@testset "native qr_full! for T = $T" for T in native_eltypes
+    rng = StableRNG(123)
+    m = 54
+    for n in (37, m, 63)
+        minmn = min(m, n)
+        A = randn(rng, T, m, n)
+        Q, R = qr_full(A, Native_HouseholderQR())
+        @test Q isa Matrix{T} && size(Q) == (m, m)
+        @test R isa Matrix{T} && size(R) == (m, n)
+        @test Q * R ≈ A
+        @test all(>=(zero(real(T))), real(diag(R)))
+        @test isunitary(Q)
+
+        Ac = similar(A)
+        Q2 = similar(Q)
+        Q2, R2 = @constinferred qr_full!(copy!(Ac, A), (Q, R), Native_HouseholderQR())
+        @test Q2 === Q
+        @test R2 === R
+        @test Q * R ≈ A
+        @test isunitary(Q)
+    end
+end

From 4dcbaf6db13c98e8bbeb7802b6767a9c86cbab44 Mon Sep 17 00:00:00 2001
From: Jutho Haegeman <jutho.haegeman@ugent.be>
Date: Tue, 4 Nov 2025 23:26:31 +0100
Subject: [PATCH 3/3] add defaults and more tests

---
 src/interface/lq.jl | 3 +++
 src/interface/qr.jl | 3 +++
 test/lq.jl          | 7 ++++++-
 test/qr.jl          | 7 ++++++-
 4 files changed, 18 insertions(+), 2 deletions(-)

diff --git a/src/interface/lq.jl b/src/interface/lq.jl
index b0d138d4..a48a4831 100644
--- a/src/interface/lq.jl
+++ b/src/interface/lq.jl
@@ -72,6 +72,9 @@ default_lq_algorithm(A; kwargs...) = default_lq_algorithm(typeof(A); kwargs...)
 function default_lq_algorithm(T::Type; kwargs...)
     throw(MethodError(default_lq_algorithm, (T,)))
 end
+function default_lq_algorithm(::Type{T}; kwargs...) where {T <: AbstractMatrix}
+    return Native_HouseholderLQ(; kwargs...)
+end
 function default_lq_algorithm(::Type{T}; kwargs...) where {T <: YALAPACK.BlasMat}
     return LAPACK_HouseholderLQ(; kwargs...)
 end
diff --git a/src/interface/qr.jl b/src/interface/qr.jl
index 7d5f8dfb..a0b7e223 100644
--- a/src/interface/qr.jl
+++ b/src/interface/qr.jl
@@ -72,6 +72,9 @@ default_qr_algorithm(A; kwargs...) = default_qr_algorithm(typeof(A); kwargs...)
 function default_qr_algorithm(T::Type; kwargs...)
     throw(MethodError(default_qr_algorithm, (T,)))
 end
+function default_qr_algorithm(::Type{T}; kwargs...) where {T <: AbstractMatrix}
+    return Native_HouseholderQR(; kwargs...)
+end
 function default_qr_algorithm(::Type{T}; kwargs...) where {T <: YALAPACK.BlasMat}
     return LAPACK_HouseholderQR(; kwargs...)
 end
diff --git a/test/lq.jl b/test/lq.jl
index 31ce2696..bd9d65b7 100644
--- a/test/lq.jl
+++ b/test/lq.jl
@@ -6,7 +6,7 @@ using LinearAlgebra: diag, I, Diagonal
 using MatrixAlgebraKit: LQViaTransposedQR, LAPACK_HouseholderQR
 
 lapack_eltypes = (Float32, Float64, ComplexF32, ComplexF64)
-native_eltypes = (lapack_eltypes..., BigFloat, Complex{BigFloat})
+native_eltypes = (lapack_eltypes..., Float16, BigFloat, Complex{BigFloat})
 
 @testset "lq_compact! for T = $T" for T in lapack_eltypes
     rng = StableRNG(123)
@@ -299,5 +299,10 @@ end
         @test Q2 === Q
         @test L * Q ≈ A
         @test isunitary(Q)
+
+        if !(T ∈ lapack_eltypes)
+            alg = MatrixAlgebraKit.select_algorithm(lq_compact!, A)
+            @test alg == Native_HouseholderLQ()
+        end
     end
 end
diff --git a/test/qr.jl b/test/qr.jl
index 6e88dfc3..ca0017ea 100644
--- a/test/qr.jl
+++ b/test/qr.jl
@@ -5,7 +5,7 @@ using StableRNGs
 using LinearAlgebra: diag, I, Diagonal
 
 lapack_eltypes = (Float32, Float64, ComplexF32, ComplexF64)
-native_eltypes = (lapack_eltypes..., BigFloat, Complex{BigFloat})
+native_eltypes = (lapack_eltypes..., Float16, BigFloat, Complex{BigFloat})
 
 @testset "qr_compact! and qr_null! for T = $T" for T in lapack_eltypes
     rng = StableRNG(123)
@@ -269,5 +269,10 @@ end
         @test R2 === R
         @test Q * R ≈ A
         @test isunitary(Q)
+
+        if !(T ∈ lapack_eltypes)
+            alg = MatrixAlgebraKit.select_algorithm(qr_compact!, A)
+            @test alg == Native_HouseholderQR()
+        end
     end
 end