Since y'all are code-literate, here the SymPy function for finding the CR-decomposition of any matrix A:
def crd(A):
"""
Computes the CR decomposition of the matrix A.
"""
rrefA, licols = A.rref() # compute RREF(A)
C = A[:, licols] # linearly indep. cols of A
r = len(licols) # = rank(A)
R = rrefA[0:r, :] # non-zero rows in RREF(A)
return C, R