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..9fa58750f 100644 --- a/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl +++ b/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...) 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..af2f799cb 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(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)) @@ -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("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 802d99bcf..36453acc0 100644 --- a/src/implementations/qr.jl +++ b/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,11 +173,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(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)) @@ -161,29 +191,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 +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) 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,