From 94103f5466a4b305766cfe0c4772e3dfe13c1e57 Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Wed, 12 Aug 2026 09:47:18 -0400
Subject: [PATCH 1/3] Extract R before constructing Q in Householder QR/LQ

cuSOLVER's `ormqr` needs a workspace of ~2*m*n elements whose size can only be
queried through a 32-bit integer, so it becomes unusable past m ~ 35000. `orgqr`
needs no workspace and is faster, but generates Q in the array holding the
reflectors, which the old order (construct Q, then copy R out of the packed A)
could not combine with returning R.

Extracting R/L first lifts the restrictions on an inplace Q (`Q === A`): it now
combines with computing R/L and with `positive = true`. The blocked path keeps
constructing Q first, as R doubles as workspace for T there. The CUSOLVER driver
now always constructs Q with `ungqr!`.

Also fixes LQ gauge fixing running with `positive = false` when L is not
computed, and makes the Native and GLA drivers reject `Q === A` instead of
silently returning garbage.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
---
 docs/src/changelog.md                         |  9 +++
 .../MatrixAlgebraKitCUDAExt.jl                |  4 ++
 ext/MatrixAlgebraKitCUDAExt/yacusolver.jl     | 22 +++++++
 ...MatrixAlgebraKitGenericLinearAlgebraExt.jl |  2 +
 src/implementations/lq.jl                     | 41 ++++++++----
 src/implementations/qr.jl                     | 64 +++++++++++++++----
 src/interface/decompositions.jl               |  4 ++
 test/decompositions/lq.jl                     |  2 +-
 test/decompositions/qr.jl                     |  2 +-
 test/testsuite/decompositions/lq.jl           | 56 ++++++++++++++++
 test/testsuite/decompositions/qr.jl           | 60 +++++++++++++++++
 11 files changed, 238 insertions(+), 28 deletions(-)

diff --git a/docs/src/changelog.md b/docs/src/changelog.md
index 58f704b6e..cc2946b68 100644
--- a/docs/src/changelog.md
+++ b/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
diff --git a/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl b/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl
index 442ee0133..a7d19dcc6 100644
--- a/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl
+++ b/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl
@@ -33,6 +33,10 @@ for f in (:geqrf!, :ungqr!, :unmqr!)
     @eval $f(::CUSOLVER, args...) = YACUSOLVER.$f(args...)
 end
 
+MatrixAlgebraKit.prefers_ungqr(::CUSOLVER) = true
+MatrixAlgebraKit.supports_unmqr(::CUSOLVER, side, trans, A, τ, C) =
+    YACUSOLVER.unmqr_worksize(side, trans, A, τ, C) > 0
+
 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...)
diff --git a/ext/MatrixAlgebraKitCUDAExt/yacusolver.jl b/ext/MatrixAlgebraKitCUDAExt/yacusolver.jl
index 0cfa64c3c..745c5274a 100644
--- a/ext/MatrixAlgebraKitCUDAExt/yacusolver.jl
+++ b/ext/MatrixAlgebraKitCUDAExt/yacusolver.jl
@@ -13,6 +13,28 @@ using CUDA.cuSOLVER: geqrf!, ormqr!, orgqr!
 const unmqr! = ormqr!
 const ungqr! = orgqr!
 
+# report a non-positive size when the 32-bit workspace query fails or overflows
+for (bname, elty) in (
+        (:cusolverDnSormqr_bufferSize, :Float32),
+        (:cusolverDnDormqr_bufferSize, :Float64),
+        (:cusolverDnCunmqr_bufferSize, :ComplexF32),
+        (:cusolverDnZunmqr_bufferSize, :ComplexF64),
+    )
+    unchecked = Symbol(:unchecked_, bname)
+    @eval function unmqr_worksize(
+            side::Char, trans::Char, A::StridedCuMatrix{$elty},
+            tau::CuVector{$elty}, C::StridedCuVecOrMat{$elty}
+        )
+        m, n = ndims(C) == 2 ? size(C) : (length(C), 1)
+        out = Ref{Cint}(0)
+        status = cuSOLVER.$unchecked(
+            cuSOLVER.dense_handle(), side, trans, m, n, length(tau),
+            A, max(1, stride(A, 2)), tau, C, max(1, stride(C, 2)), out
+        )
+        return status == cuSOLVER.CUSOLVER_STATUS_SUCCESS ? Int(out[]) : -1
+    end
+end
+
 # Wrapper for SVD via QR Iteration
 for (bname, fname, elty, relty) in
     (
diff --git a/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl b/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl
index 88d043350..93f846c81 100644
--- a/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl
+++ b/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)
diff --git a/src/implementations/lq.jl b/src/implementations/lq.jl
index bef88d4f8..dd387c1e3 100644
--- a/src/implementations/lq.jl
+++ b/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,9 +151,16 @@ 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("inplace Q only supported if matrix is wide (`m <= n`) and using the unblocked algorithm (`blocksize = 1`)"))
+        (computeL && Base.mightalias(L, A)) &&
+            throw(ArgumentError("inplace Q 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))
@@ -154,24 +168,22 @@ function lq_householder!(
             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("inplace Q 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')
diff --git a/src/implementations/qr.jl b/src/implementations/qr.jl
index 802d99bcf..637919aee 100644
--- a/src/implementations/qr.jl
+++ b/src/implementations/qr.jl
@@ -125,6 +125,34 @@ for f in (:geqrt!, :gemqrt!, :geqp3!, :geqrf!, :ungqr!, :unmqr!)
     end
 end
 
+# cuSOLVER's 32-bit ormqr workspace query breaks down for large matrices
+supports_unmqr(::Driver, side, trans, A, τ, C) = true
+
+# cuSOLVER generates Q faster with ungqr! than by applying the reflectors with unmqr!
+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)
+    (prefers_ungqr(driver) || !supports_unmqr(driver, 'L', 'N', A, τ, Q)) &&
+        return _qr_buildQ_ungqr!(driver, Q, A, τ, minmn)
+    return unmqr!(driver, 'L', 'N', A, τ, one!(Q))
+end
+
+# build Q in its own space: copy in the reflectors, unit vectors elsewhere
+function _qr_buildQ_ungqr!(driver::Driver, Q::AbstractMatrix, A::AbstractMatrix, τ, minmn::Int)
+    size(Q, 2) > minmn && one!(Q)
+    copyto!(view(Q, :, 1:minmn), view(A, :, 1:minmn))
+    ungqr!(driver, Q, τ)
+    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,11 +177,17 @@ 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("inplace Q only supported if matrix is tall (`m >= n`) and using the unblocked algorithm (`blocksize = 1`)"))
+        (computeR && Base.mightalias(R, A)) &&
+            throw(ArgumentError("inplace Q 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))
@@ -161,29 +195,27 @@ function qr_householder!(
             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 +227,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"))
     # positive = true regardless of setting
 
     m, n = size(A)
@@ -257,6 +291,8 @@ function qr_null_householder!(
         N = gemqrt!(driver, 'L', 'N', A, T, N)
     else
         A, τ = geqrf!(driver, A)
+        supports_unmqr(driver, 'L', 'N', A, τ, N) ||
+            throw(ArgumentError(lazy"$driver cannot construct the null space for these dimensions"))
         N = unmqr!(driver, 'L', 'N', A, τ, N)
     end
     return N
diff --git a/src/interface/decompositions.jl b/src/interface/decompositions.jl
index 4077214e8..1d423564a 100644
--- a/src/interface/decompositions.jl
+++ b/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(;
diff --git a/test/decompositions/lq.jl b/test/decompositions/lq.jl
index c5c2f047a..2b13ee195 100644
--- a/test/decompositions/lq.jl
+++ b/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
diff --git a/test/decompositions/qr.jl b/test/decompositions/qr.jl
index 666d641c5..7601237fc 100644
--- a/test/decompositions/qr.jl
+++ b/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
diff --git a/test/testsuite/decompositions/lq.jl b/test/testsuite/decompositions/lq.jl
index d53ff4438..d50691bbb 100644
--- a/test/testsuite/decompositions/lq.jl
+++ b/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),
diff --git a/test/testsuite/decompositions/qr.jl b/test/testsuite/decompositions/qr.jl
index d5ebf402f..a8dbd9376 100644
--- a/test/testsuite/decompositions/qr.jl
+++ b/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,

From e16420079efc2f0ff28b6c61812c71d732f0feca Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Wed, 12 Aug 2026 10:09:52 -0400
Subject: [PATCH 2/3] Drop the cuSOLVER workspace queries

They are no longer used now that the CUSOLVER driver always constructs Q with
`ungqr!`. Inline the ungqr!-based Q construction into `_qr_buildQ!`.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
---
 .../MatrixAlgebraKitCUDAExt.jl                |  2 --
 ext/MatrixAlgebraKitCUDAExt/yacusolver.jl     | 22 ----------------
 src/implementations/qr.jl                     | 26 +++++++------------
 3 files changed, 10 insertions(+), 40 deletions(-)

diff --git a/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl b/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl
index a7d19dcc6..9fa58750f 100644
--- a/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl
+++ b/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl
@@ -34,8 +34,6 @@ for f in (:geqrf!, :ungqr!, :unmqr!)
 end
 
 MatrixAlgebraKit.prefers_ungqr(::CUSOLVER) = true
-MatrixAlgebraKit.supports_unmqr(::CUSOLVER, side, trans, A, τ, C) =
-    YACUSOLVER.unmqr_worksize(side, trans, A, τ, C) > 0
 
 MatrixAlgebraKit.supports_svd_full(::CUSOLVER, f::Symbol) = f in (:qr_iteration, :jacobi, :svd_polar)
 
diff --git a/ext/MatrixAlgebraKitCUDAExt/yacusolver.jl b/ext/MatrixAlgebraKitCUDAExt/yacusolver.jl
index 745c5274a..0cfa64c3c 100644
--- a/ext/MatrixAlgebraKitCUDAExt/yacusolver.jl
+++ b/ext/MatrixAlgebraKitCUDAExt/yacusolver.jl
@@ -13,28 +13,6 @@ using CUDA.cuSOLVER: geqrf!, ormqr!, orgqr!
 const unmqr! = ormqr!
 const ungqr! = orgqr!
 
-# report a non-positive size when the 32-bit workspace query fails or overflows
-for (bname, elty) in (
-        (:cusolverDnSormqr_bufferSize, :Float32),
-        (:cusolverDnDormqr_bufferSize, :Float64),
-        (:cusolverDnCunmqr_bufferSize, :ComplexF32),
-        (:cusolverDnZunmqr_bufferSize, :ComplexF64),
-    )
-    unchecked = Symbol(:unchecked_, bname)
-    @eval function unmqr_worksize(
-            side::Char, trans::Char, A::StridedCuMatrix{$elty},
-            tau::CuVector{$elty}, C::StridedCuVecOrMat{$elty}
-        )
-        m, n = ndims(C) == 2 ? size(C) : (length(C), 1)
-        out = Ref{Cint}(0)
-        status = cuSOLVER.$unchecked(
-            cuSOLVER.dense_handle(), side, trans, m, n, length(tau),
-            A, max(1, stride(A, 2)), tau, C, max(1, stride(C, 2)), out
-        )
-        return status == cuSOLVER.CUSOLVER_STATUS_SUCCESS ? Int(out[]) : -1
-    end
-end
-
 # Wrapper for SVD via QR Iteration
 for (bname, fname, elty, relty) in
     (
diff --git a/src/implementations/qr.jl b/src/implementations/qr.jl
index 637919aee..b2f0410fc 100644
--- a/src/implementations/qr.jl
+++ b/src/implementations/qr.jl
@@ -125,10 +125,8 @@ for f in (:geqrt!, :gemqrt!, :geqp3!, :geqrf!, :ungqr!, :unmqr!)
     end
 end
 
-# cuSOLVER's 32-bit ormqr workspace query breaks down for large matrices
-supports_unmqr(::Driver, side, trans, A, τ, C) = true
-
-# cuSOLVER generates Q faster with ungqr! than by applying the reflectors with unmqr!
+# 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
@@ -140,16 +138,14 @@ function _qr_copyR!(R::AbstractMatrix, A::AbstractMatrix, jpvt = nothing)
 end
 
 function _qr_buildQ!(driver::Driver, Q::AbstractMatrix, A::AbstractMatrix, τ, minmn::Int)
-    (prefers_ungqr(driver) || !supports_unmqr(driver, 'L', 'N', A, τ, Q)) &&
-        return _qr_buildQ_ungqr!(driver, Q, A, τ, minmn)
-    return unmqr!(driver, 'L', 'N', A, τ, one!(Q))
-end
-
-# build Q in its own space: copy in the reflectors, unit vectors elsewhere
-function _qr_buildQ_ungqr!(driver::Driver, Q::AbstractMatrix, A::AbstractMatrix, τ, minmn::Int)
-    size(Q, 2) > minmn && one!(Q)
-    copyto!(view(Q, :, 1:minmn), view(A, :, 1:minmn))
-    ungqr!(driver, Q, τ)
+    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
 
@@ -291,8 +287,6 @@ function qr_null_householder!(
         N = gemqrt!(driver, 'L', 'N', A, T, N)
     else
         A, τ = geqrf!(driver, A)
-        supports_unmqr(driver, 'L', 'N', A, τ, N) ||
-            throw(ArgumentError(lazy"$driver cannot construct the null space for these dimensions"))
         N = unmqr!(driver, 'L', 'N', A, τ, N)
     end
     return N

From 172e4a3eb11ea7077404de9cea49395e44d26294 Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Thu, 13 Aug 2026 03:45:39 -0400
Subject: [PATCH 3/3] code review

---
 src/implementations/lq.jl | 6 +++---
 src/implementations/qr.jl | 6 +++---
 2 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/src/implementations/lq.jl b/src/implementations/lq.jl
index dd387c1e3..af2f799cb 100644
--- a/src/implementations/lq.jl
+++ b/src/implementations/lq.jl
@@ -154,9 +154,9 @@ function lq_householder!(
     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("inplace Q only supported if matrix is wide (`m <= n`) and using the unblocked algorithm (`blocksize = 1`)"))
+            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("inplace Q only supported if L does not share memory with A"))
+            throw(ArgumentError("in-place Q is only supported if L does not share memory with A"))
     end
 
     if blocksize > 1
@@ -319,7 +319,7 @@ function lq_via_qr!(
     )
     # 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("inplace Q only supported if L does not share memory with 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')
diff --git a/src/implementations/qr.jl b/src/implementations/qr.jl
index b2f0410fc..36453acc0 100644
--- a/src/implementations/qr.jl
+++ b/src/implementations/qr.jl
@@ -176,9 +176,9 @@ function qr_householder!(
     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("inplace Q only supported if matrix is tall (`m >= n`) and using the unblocked algorithm (`blocksize = 1`)"))
+            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("inplace Q only supported if R does not share memory with A"))
+            throw(ArgumentError("in-place Q is only supported if R does not share memory with A"))
     end
 
     # Compute QR in packed form
@@ -224,7 +224,7 @@ function qr_householder!(
     pivoted &&
         throw(ArgumentError(lazy"$driver does not provide a pivoted QR decomposition"))
     Q === A &&
-        throw(ArgumentError(lazy"$driver does not provide an inplace Q"))
+        throw(ArgumentError(lazy"$driver does not provide an in-place Q"))
     # positive = true regardless of setting
 
     m, n = size(A)