using LinearAlgebra
using SparseArrays

function R(n,p,q,c,s)
    M=sparse(Float64,I,n,n)
    M[p,p]=c; M[p,q]=s
    M[q,p]=-s; M[q,q]=c
    return M
end

function givens(A,b)
    m,n=size(A)
    if m<n
        error("wide matrices not supported!\n")
    end
    for p=1:n
        for q=m:-1:p+1
            r1=A[p,p]; r2=A[q,p]
            z=sqrt(r1*r1+r2*r2)
            if z>0
                c=r1/z
                s=r2/z
            else
                c=1
                s=0
            end
            M=R(m,p,q,c,s)
            A=M*A
            b=M*b
        end
    end
    return A,b
end