Convert between QR and QRCompactWY

[sethaxen]

This is certainly niche, but I have a use case where I manually constructed a QR object but then wanted a QRCompactWY object to use it most efficiently. It's possible to interconvert between the representations used by the 2 objects, but no convert methods exist for this. Would the following implementations be useful to add to LinearAlgebra?

function Base.convert(::Type{S}, F::LinearAlgebra.QR; blocksize::Int=36) where {S<:LinearAlgebra.QRCompactWY}
    factors = Base.copymutable(F.factors)
    k = length(F.τ)
    T = similar(factors, min(blocksize, k), k)
    _generate_compact_WY!(factors, F.τ, T)
    return S(factors, T)
end

function _generate_compact_WY!(V::AbstractMatrix, τ::AbstractVector, T::AbstractMatrix)
    m, _ = size(V)
    k = length(τ)
    nb = size(T, 1)
    for i in 1:nb:k
        nbi = min(k - i + 1, nb)
        # select same blocks as LAPACK xGEQRT
        @views _generate_compact_WY_block!(V[i:m, i:i+nbi-1], τ[i:i+nbi-1], T[1:nbi, i:i+nbi-1])
    end
    return T
end

# fill block of T in the same way as LAPACK xGEQRT2
function _generate_compact_WY_block!(V::AbstractMatrix, τ::AbstractVector, T::AbstractMatrix)
    m, n = size(V)
    T[1, 1] = τ[1]
    for i in 2:n
        T[i,i] = τ[i]
        aii = V[i,i]
        V[i,i] = 1
        ti = view(T, 1:i-1, i)
        mul!(ti, view(V, i:m, 1:i-1)', view(V, i:m, i), -τ[i], 0)
        V[i,i] = aii
        lmul!(UpperTriangular(view(T, 1:i-1, 1:i-1)), ti)
    end
    return T
end

function Base.convert(::Type{S}, F::LinearAlgebra.QRCompactWY) where {S<:LinearAlgebra.QR}
    nb, k = size(F.T)
    τ = similar(F.T, (axes(F.T, 2),))
    for i in 1:nb:k
        nbi = min(k - i + 1, nb)
        for (ib, j) in enumerate(i:i+nbi-1)
            τ[j] = F.T[ib, j]
        end
    end
    return S(copy(F.factors), τ)
end

Here's a quick demo:

julia> A = randn(40, 100); Fcomp = qr(A); F = LinearAlgebra.qrfactUnblocked!(copy(A)); 

julia> convert(LinearAlgebra.QR, Fcomp).τ ≈ F.τ
true

julia> convert(LinearAlgebra.QR, convert(LinearAlgebra.QRCompactWY, F)).τ == F.τ
true

julia> Matrix(convert(LinearAlgebra.QR, Fcomp)) ≈ Matrix(convert(LinearAlgebra.QRCompactWY, F)) ≈ A
true