QR decomposition by Givens rotations

qr-givens.jl

"""
QR decomposition by givens rotations of matrix without pivoting.

```math
A = Q R
```

Thin (reduced) method will produce `Q` and `R` in truncated form,
in opposite case `thin=false` Q is full, but R is still reduced, see [`qr`](@ref).
"""
@generated function qr_givens_unrolled(::Size{sA}, A::StaticMatrix{<:Any, <:Any, TA}, thin::Type) where {sA, TA}
    mQ = nQ = mR = m = sA[1]
    nR = n = sA[2]
    m > n && (mR = n)
    m > n && thin <: Type{Val{true}} && (nQ = n)

    Q = [Symbol("Q_$(i)_$(j)") for i = 1:m, j = 1:m]
    R = [Symbol("R_$(i)_$(j)") for i = 1:m, j = 1:n]

    initQ = [:($(Q[i, j]) = $(i == j ? one : zero)(T)) for i = 1:m, j = 1:m]  # Q .= eye(A)
    initR = [:($(R[i, j]) = T(A[$i, $j])) for i = 1:m, j = 1:n]               # R .= A

    #for j = 1:n-1
    #    for i = m:-1:j+1
    #        (c, s, r) = LinAlg.givensAlgorithm(R[i-1,j], R[i,j])
    #        # R = G'*R, Q = Q*G
    #        R[i-1,j] = r
    #        R[i  ,j] = zero(T)
    #        for k = j+1:n
    #            r1       =  R[i-1,k]
    #            r2       =  R[i  ,k]
    #            R[i-1,k] =  c*r1 + s*r2
    #            R[i  ,k] = -s*r1 + c*r2
    #        end
    #        for k = 1:m
    #            r1       =  Q[k,i-1]
    #            r2       =  Q[k,i]
    #            Q[k,i-1] =  c*r1 + s*r2
    #            Q[k,i]   = -s*r1 + c*r2
    #        end
    #    end
    #end
    code = quote end
    for j = 1:n
        for i = m:-1:j+1
            push!(code.args, :((c, s, r) = LinAlg.givensAlgorithm($(R[i-1,j]), $(R[i,j]))))
            push!(code.args, :($(R[i-1,j]) =  r))
            push!(code.args, :($(R[i  ,j]) =  zero(T)))
            for k = j+1:n
                push!(code.args, :(r1          =  $(R[i-1,k])))
                push!(code.args, :(r2          =  $(R[i  ,k])))
                push!(code.args, :($(R[i-1,k]) =  c*r1 + s*r2))
                push!(code.args, :($(R[i  ,k]) = -s*r1 + c*r2))
            end
            for k = 1:m
                push!(code.args, :(r1          =  $(Q[k,i-1])))
                push!(code.args, :(r2          =  $(Q[k,i])))
                push!(code.args, :($(Q[k,i-1]) =  c*r1 + s*r2))
                push!(code.args, :($(Q[k,i])   = -s*r1 + c*r2))
            end
        end
    end

    return quote
        @_inline_meta
        T = promote_op(matprod, promote_type(typeof(sqrt(one(TA))),Float32), arithmetic_closure(TA))
        @inbounds $(Expr(:block, initQ...))
        @inbounds $(Expr(:block, initR...))
        @inbounds $code
        @inbounds return similar_type(A, T, $(Size(mQ,nQ)))(tuple($(Q[1:mQ,1:nQ]...))),
                         similar_type(A, T, $(Size(mR,nR)))(tuple($(R[1:mR,1:nR]...)))
    end
end