K = RR

# Input matrix here
A = matrix(K, [[-3,     32/5,          4   ],
               [4,      24/5,          3   ],
               [5,   6*sqrt(2),  5*sqrt(2)]])

m,n = (3,3)


# Calculate r, c, s
def givens(a,b):
    if b == 0:
        c = 1
        s = 0
    else:
        if abs(b) > abs(a):
            r = a / b;
            s = 1 / sqrt(1 + r**2)
            c = s*r
        else:
            r = b / a
            c = 1 / sqrt(1 + r**2)
            s = c*r
    return r, c, s

Q = matrix(K,m,n,1)
R = copy(A)

# Calculate QR
for j in range(0,n):
    for i in range(m-1, j, -1):
        r, c, s = givens(R[i-1,j], R[i,j])
        R[i-1:i+1,j:n] = matrix([[c,s],[-s,c]]) * R[i-1:i+1,j:n]
        Q[i-1:i+1,:] = matrix([[c,s],[-s,c]]) * Q[i-1:i+1,:]


Q = Q.transpose()

print "-" * 10
print Q
print
print R
print
print Q * R