Extract R before constructing Q in Householder QR/LQ

docs/src/changelog.md

@@ -24,12 +24,21 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
 
 ### Changed
 
+- `qr_compact!`, `qr_full!`, `lq_compact!` and `lq_full!` now extract `R` (or `L`) before
+  constructing `Q`, so that an inplace `Q` (supplying `A` itself as output for `Q`) can be combined
+  with computing `R` (or `L`) and with `positive = true`.
+- The CUSOLVER driver constructs `Q` with `ungqr!` instead of `unmqr!`, which is both faster and
+  avoids the large workspace of `ormqr`, whose 32-bit size query fails altogether for large
+  matrices.
+
 ### Deprecated
 
 ### Removed
 
 ### Fixed
 
+- LQ decompositions no longer gauge fix `Q` when `positive = false` and `L` is not computed.
+
 ### Performance
 
 ## [0.6.9](https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/compare/v0.6.8...v0.6.9) - 2026-07-10

ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl

@@ -33,6 +33,8 @@ for f in (:geqrf!, :ungqr!, :unmqr!)
     @eval $f(::CUSOLVER, args...) = YACUSOLVER.$f(args...)
 end
 
+MatrixAlgebraKit.prefers_ungqr(::CUSOLVER) = true
+
 MatrixAlgebraKit.supports_svd_full(::CUSOLVER, f::Symbol) = f in (:qr_iteration, :jacobi, :svd_polar)
 
 function gesvd!(::CUSOLVER, A::StridedCuMatrix, S::StridedCuVector, U::StridedCuMatrix, Vᴴ::StridedCuMatrix; kwargs...)

ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl

@@ -61,6 +61,8 @@ function qr_householder!(
         throw(ArgumentError(lazy"$driver does not provide a blocked QR decomposition"))
     pivoted &&
         throw(ArgumentError(lazy"$driver does not provide a pivoted QR decomposition"))
+    Q === A &&
+        throw(ArgumentError(lazy"$driver does not provide an inplace Q"))
 
     m, n = size(A)
     minmn = min(m, n)

src/implementations/lq.jl

@@ -123,6 +123,13 @@ for f in (:gelqt!, :gemlqt!, :gelqf!, :unglq!, :unmlq!)
     end
 end
 
+# copy L out of the packed factorization, leaving the reflectors in A intact
+function _lq_copyL!(L::AbstractMatrix, A::AbstractMatrix)
+    copyto!(L, view(A, axes(L)...))
+    lowertriangular!(L)
+    return L
+end
+
 @inline lq_householder!(A, L, Q; driver::Driver = DefaultDriver(), kwargs...) =
     lq_householder!(driver, A, L, Q; kwargs...)
 lq_householder!(::DefaultDriver, A, L, Q; kwargs...) =
@@ -144,34 +151,39 @@ function lq_householder!(
     computeL = length(L) > 0
     inplaceQ = Q === A
 
-    (inplaceQ && (computeL || positive || blocksize > 1 || n < m)) &&
-        throw(ArgumentError("inplace Q only supported if matrix is wide (`m <= n`), L is not required, and using the unblocked algorithm (`blocksize = 1`) with `positive = false`"))
+    if inplaceQ
+        # unglq! builds Q in the space of A, so L has to be extracted first and cannot alias A
+        (blocksize == 1 && n >= m) ||
+            throw(ArgumentError(lazy"in-place Q is only supported if matrix is wide (`$m ≤ $n`) and using the unblocked algorithm (`blocksize = $blocksize`)"))
+        (computeL && Base.mightalias(L, A)) &&
+            throw(ArgumentError("in-place Q is only supported if L does not share memory with A"))
+    end
+
     if blocksize > 1
+        # L doubles as workspace for T, so Q is constructed before L is extracted
         mb = min(minmn, blocksize)
         if computeL # first use L as space for T
             A, T = gelqt!(driver, A, view(L, 1:mb, 1:minmn))
         else
             A, T = gelqt!(driver, A, similar(A, mb, minmn))
         end
         Q = gemlqt!(driver, 'R', 'N', A, T, one!(Q))
+        computeL && _lq_copyL!(L, A)
+        positive && gaugefix!(lq_householder!, computeL ? L : nothing, Q, diagview(A))
     else
         A, τ = gelqf!(driver, A)
+        computeL && _lq_copyL!(L, A)
+        Lf = computeL ? L : nothing
         if inplaceQ
-            Q = unglq!(driver, A, τ)
+            Ld = positive ? copy(diagview(A)) : nothing # unglq! destroys the diagonal of A
+            unglq!(driver, A, τ) # Q === A, so no need to rebind Q
+            positive && gaugefix!(lq_householder!, Lf, Q, Ld)
         else
             Q = unmlq!(driver, 'R', 'N', A, τ, one!(Q))
+            positive && gaugefix!(lq_householder!, Lf, Q, diagview(A))
         end
     end
 
-    if computeL
-        # we need to first copy then gaugefix - avoiding aliasing between L and Ld for broadcast
-        Ld = diagview(A)
-        copyto!(L, lowertriangular!(view(A, axes(L)...)))
-        positive && gaugefix!(lq_householder!, L, Q, Ld)
-    else
-        gaugefix!(lq_householder!, nothing, Q, diagview(A))
-    end
-
     return L, Q
 end
 function lq_householder!(
@@ -183,6 +195,8 @@ function lq_householder!(
         throw(ArgumentError(lazy"$driver does not provide a blocked LQ decomposition"))
     pivoted &&
         throw(ArgumentError(lazy"$driver does not provide a pivoted LQ decomposition"))
+    Q === A &&
+        throw(ArgumentError(lazy"$driver does not provide an inplace Q"))
     # positive = true regardless of setting
 
     m, n = size(A)
@@ -303,6 +317,9 @@ end
 function lq_via_qr!(
         A::AbstractMatrix, L::AbstractMatrix, Q::AbstractMatrix, qr_alg::AbstractAlgorithm
     )
+    # Q is written before L, so an L that aliases A would corrupt an inplace Q
+    (A === Q && !isempty(L) && Base.mightalias(L, A)) &&
+        throw(ArgumentError("in-place Q is only supported if L does not share memory with A"))
     At = adjoint!(similar(A'), A)::AbstractMatrix
     Qt = (A === Q) ? At : similar(Q')
     Lt = similar(L')

src/implementations/qr.jl

@@ -125,6 +125,30 @@ for f in (:geqrt!, :gemqrt!, :geqp3!, :geqrf!, :ungqr!, :unmqr!)
     end
 end
 
+# cuSOLVER generates Q faster with ungqr! than by applying the reflectors with unmqr!,
+# and avoids the large workspace of ormqr, whose 32-bit size query fails for large matrices
+prefers_ungqr(::Driver) = false
+
+# copy R out of the packed factorization, leaving the reflectors in A intact
+function _qr_copyR!(R::AbstractMatrix, A::AbstractMatrix, jpvt = nothing)
+    Rp = isnothing(jpvt) ? R : view(R, :, jpvt)
+    copyto!(Rp, view(A, axes(R)...))
+    uppertriangular!(Rp)
+    return R
+end
+
+function _qr_buildQ!(driver::Driver, Q::AbstractMatrix, A::AbstractMatrix, τ, minmn::Int)
+    if prefers_ungqr(driver)
+        # build Q in its own space: copy in the reflectors, unit vectors elsewhere
+        size(Q, 2) > minmn && one!(Q)
+        copyto!(view(Q, :, 1:minmn), view(A, :, 1:minmn))
+        ungqr!(driver, Q, τ)
+    else
+        Q = unmqr!(driver, 'L', 'N', A, τ, one!(Q))
+    end
+    return Q
+end
+
 @inline qr_householder!(A, Q, R; driver::Driver = DefaultDriver(), kwargs...) =
     qr_householder!(driver, A, Q, R; kwargs...)
 qr_householder!(::DefaultDriver, A, Q, R; kwargs...) =
@@ -149,41 +173,45 @@ function qr_householder!(
     computeR = length(R) > 0
     inplaceQ = Q === A
 
-    (inplaceQ && (computeR || positive || blocksize > 1 || m < n)) &&
-        throw(ArgumentError("inplace Q only supported if matrix is tall (`m >= n`), R is not required, and using the unblocked algorithm (`blocksize = 1`) with `positive = false`"))
+    if inplaceQ
+        # ungqr! builds Q in the space of A, so R has to be extracted first and cannot alias A
+        (blocksize == 1 && m >= n) ||
+            throw(ArgumentError(lazy"in-place Q is only supported if matrix is tall (`$m >= $n`) and using the unblocked algorithm (`blocksize = $blocksize`)"))
+        (computeR && Base.mightalias(R, A)) &&
+            throw(ArgumentError("in-place Q is only supported if R does not share memory with A"))
+    end
 
     # Compute QR in packed form
     if blocksize > 1
+        # R doubles as workspace for T, so Q is constructed before R is extracted
         nb = min(minmn, blocksize)
         if computeR # first use R as space for T
             A, T = geqrt!(driver, A, view(R, 1:nb, 1:minmn))
         else
             A, T = geqrt!(driver, A, similar(A, nb, minmn))
         end
         Q = gemqrt!(driver, 'L', 'N', A, T, one!(Q))
+        computeR && _qr_copyR!(R, A)
+        positive && gaugefix!(qr_householder!, Q, computeR ? R : nothing, diagview(A))
     else
         if pivoted
             A, τ, jpvt = geqp3!(driver, A)
+            computeR && _qr_copyR!(R, A, jpvt)
         else
             A, τ = geqrf!(driver, A)
+            computeR && _qr_copyR!(R, A)
         end
+        Rf = computeR ? R : nothing # gaugefix! rescales rows, which commutes with the pivoting
         if inplaceQ
-            Q = ungqr!(driver, A, τ)
+            Rd = positive ? copy(diagview(A)) : nothing # ungqr! destroys the diagonal of A
+            ungqr!(driver, A, τ) # Q === A, so no need to rebind Q
+            positive && gaugefix!(qr_householder!, Q, Rf, Rd)
         else
-            Q = unmqr!(driver, 'L', 'N', A, τ, one!(Q))
+            Q = _qr_buildQ!(driver, Q, A, τ, minmn)
+            positive && gaugefix!(qr_householder!, Q, Rf, diagview(A))
         end
     end
 
-    if computeR
-        # we need to first copy then gaugefix - avoiding aliasing between R and Rd for broadcast
-        Rd = diagview(A)
-        Rf = pivoted ? view(R, :, jpvt) : R
-        copyto!(Rf, uppertriangular!(view(A, axes(R)...)))
-        positive && gaugefix!(qr_householder!, Q, Rf, Rd)
-    elseif positive
-        gaugefix!(qr_householder!, Q, nothing, diagview(A))
-    end
-
     return Q, R
 end
 function qr_householder!(
@@ -195,6 +223,8 @@ function qr_householder!(
         throw(ArgumentError(lazy"$driver does not provide a blocked QR decomposition"))
     pivoted &&
         throw(ArgumentError(lazy"$driver does not provide a pivoted QR decomposition"))
+    Q === A &&
+        throw(ArgumentError(lazy"$driver does not provide an in-place Q"))
     # positive = true regardless of setting
 
     m, n = size(A)

src/interface/decompositions.jl

@@ -76,6 +76,10 @@ Algorithm type to denote the algorithm for computing QR, RQ, QL or LQ decomposit
 
 Depending on the driver, various other keywords may be (un)available to customize the implementation.
 The optional `driver` keyword can be used to choose between different implementations of this algorithm.
+
+`Q` may be computed in the space of the input matrix `A`, by supplying `A` itself as output for
+`Q`. This requires `blocksize = 1`, a tall `A` for QR (a wide `A` for LQ), and `R` (or `L`) to not
+share memory with `A`.
 """
 @algdef Householder
 function Householder(;

test/decompositions/lq.jl

@@ -57,7 +57,7 @@ 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)
+            TestSuite.test_lq(T, (m, n); test_pivoted = false, test_blocksize = false, test_inplaceQ = false)
             GENERIC_LQ_ALGS = (Householder(; driver = Native()), LQViaTransposedQR(Householder(; driver = GLA())))
             TestSuite.test_lq_algs(T, (m, n), GENERIC_LQ_ALGS)
         end

test/decompositions/qr.jl

@@ -56,7 +56,7 @@ 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)
+            TestSuite.test_qr(T, (m, n); test_pivoted = false, test_blocksize = false, test_inplaceQ = false)
             GENERIC_QR_ALGS = (Householder(; driver = Native()), Householder(; driver = GLA()))
             TestSuite.test_qr_algs(T, (m, n), GENERIC_QR_ALGS)
         end

test/testsuite/decompositions/lq.jl

@@ -6,6 +6,7 @@ function test_lq(T::Type, sz; test_null = true, kwargs...)
         test_lq_compact(T, sz; kwargs...)
         test_lq_full(T, sz; kwargs...)
         test_null && test_lq_null(T, sz; kwargs...)
+        test_lq_inplaceQ(T, sz; kwargs...)
     end
 end
 
@@ -75,6 +76,61 @@ function test_lq_compact(
     end
 end
 
+# test using `A` itself as output for `Q`
+function test_lq_inplaceQ(
+        T::Type, sz;
+        test_inplaceQ = true, test_positive = true,
+        atol::Real = 0, rtol::Real = precision(T),
+        kwargs...
+    )
+    (sz isa Tuple && length(sz) == 2) || return nothing
+    m, n = sz
+    n >= m || return nothing # inplace Q requires a wide matrix
+    summary_str = testargs_summary(T, sz)
+    return @testset "lq_compact! inplace Q $summary_str" begin
+        A = instantiate_matrix(T, sz)
+        L = similar(A, (m, m))
+
+        if !test_inplaceQ
+            Ain = deepcopy(A)
+            @test_throws Exception lq_compact!(Ain, (L, Ain))
+        else
+            for positive in (test_positive ? (false, true) : (false,))
+                Ain = deepcopy(A)
+                L2, Q = lq_compact!(Ain, (L, Ain); positive)
+                @test Q === Ain
+                @test L2 * Q ≈ A
+                @test isisometric(Q; side = :right, atol, rtol)
+                @test istril(L2)
+                if positive
+                    @test has_positive_diagonal(L2)
+                end
+
+                # L is not required
+                Ain2 = deepcopy(A)
+                _, Q2 = lq_compact!(Ain2, (similar(A, (0, 0)), Ain2); positive)
+                @test Q2 === Ain2
+                @test Q2 ≈ Q
+            end
+
+            if m == n
+                Ain = deepcopy(A)
+                Lf, Q = lq_full!(Ain, (similar(A, (m, n)), Ain))
+                @test Q === Ain
+                @test Lf * Q ≈ A
+                @test isunitary(Q; atol, rtol)
+                @test has_positive_diagonal(Lf)
+            end
+
+            # blocked algorithm and aliased L are not supported
+            Ain = deepcopy(A)
+            @test_throws ArgumentError lq_compact!(Ain, (L, Ain); blocksize = 2)
+            Ain = deepcopy(A)
+            @test_throws ArgumentError lq_compact!(Ain, (view(Ain, :, 1:m), Ain))
+        end
+    end
+end
+
 function test_lq_compact_algs(
         T::Type, sz, algs;
         atol::Real = 0, rtol::Real = precision(T),

test/testsuite/decompositions/qr.jl

@@ -6,6 +6,7 @@ function test_qr(T::Type, sz; kwargs...)
         test_qr_compact(T, sz; kwargs...)
         test_qr_full(T, sz; kwargs...)
         test_qr_null(T, sz; kwargs...)
+        test_qr_inplaceQ(T, sz; kwargs...)
     end
 end
 
@@ -112,6 +113,65 @@ function test_qr_compact_algs(
     end
 end
 
+# test using `A` itself as output for `Q`
+function test_qr_inplaceQ(
+        T::Type, sz;
+        test_inplaceQ = true, test_positive = true, test_pivoted = true,
+        atol::Real = 0, rtol::Real = precision(T),
+        kwargs...
+    )
+    (sz isa Tuple && length(sz) == 2) || return nothing
+    m, n = sz
+    m >= n || return nothing # inplace Q requires a tall matrix
+    summary_str = testargs_summary(T, sz)
+    return @testset "qr_compact! inplace Q $summary_str" begin
+        A = instantiate_matrix(T, sz)
+        R = similar(A, (n, n))
+
+        if !test_inplaceQ
+            Ain = deepcopy(A)
+            @test_throws Exception qr_compact!(Ain, (Ain, R))
+        else
+            for positive in (test_positive ? (false, true) : (false,)),
+                    pivoted in (test_pivoted ? (false, true) : (false,))
+
+                Ain = deepcopy(A)
+                Q, R2 = qr_compact!(Ain, (Ain, R); positive, pivoted)
+                @test Q === Ain
+                @test Q * R2 ≈ A
+                @test isisometric(Q; atol, rtol)
+                if !pivoted
+                    @test istriu(R2)
+                    if positive
+                        @test has_positive_diagonal(R2)
+                    end
+                end
+
+                # R is not required
+                Ain2 = deepcopy(A)
+                Q2, = qr_compact!(Ain2, (Ain2, similar(A, (0, 0))); positive, pivoted)
+                @test Q2 === Ain2
+                @test Q2 ≈ Q
+            end
+
+            if m == n
+                Ain = deepcopy(A)
+                Q, Rf = qr_full!(Ain, (Ain, similar(A, (m, n))))
+                @test Q === Ain
+                @test Q * Rf ≈ A
+                @test isunitary(Q; atol, rtol)
+                @test has_positive_diagonal(Rf)
+            end
+
+            # blocked algorithm and aliased R are not supported
+            Ain = deepcopy(A)
+            @test_throws ArgumentError qr_compact!(Ain, (Ain, R); blocksize = 2)
+            Ain = deepcopy(A)
+            @test_throws ArgumentError qr_compact!(Ain, (Ain, view(Ain, 1:n, :)))
+        end
+    end
+end
+
 function test_qr_full(
         T::Type, sz;
         test_positive = true, test_pivoted = true, test_blocksize = true,