Skip to content

Instantly share code, notes, and snippets.

@denius

denius/qr-givens.jl

Created Oct 29, 2017
Embed
What would you like to do?
QR decomposition by Givens rotations
"""
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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.