def back_substitution(M, X, n):
    """Compute MR = X for some n*n upper-triangular matrix M and
       some n*m matrix X."""
    R = [[0] * len(X[0]) for _ in xrange(0, n)]
        
    def substitute(equation, S, row, col):
        """Substitute a solved value into all equations above it
           in the matrix."""
        for r in xrange(0, row):
            S[r] += M[r][col] * R[col][equation]
            
    def solve(equation, S, row, col):
        """Solve an equation of the from ax + b = c, where a=M[row][col],
           b=S[col], and c = X[row][equation]."""
        if 0 != M[row][col]:
            R[col][equation] = (X[row][equation] - S[col]) / M[row][col]
        else:
            R[col][equation] = 0
                
        if row > 0 and col > 0:
            substitute(equation, S, row, col)
            solve(equation, S, row - 1, col - 1)
    
    # perform back-substitution on each column vector
    for i in xrange(0, len(X[0])):
        S = [0] * n
        solve(i, S, n - 1, n - 1)
    
    return R

M = [
    [1, -2, 1],
    [0, 1, 6],
    [0, 0, 1]
]

M_inv = back_substitution(
    M, 
    [
        [1, 0, 0], 
        [0, 1, 0], 
        [0, 0, 1]
    ], 
    3
)

def mult(A, B):
    """Multiply two n*n matrices."""
    rng = xrange(0, len(A))
    R = [[0] * len(A) for _ in rng]
    for i in rng:
        for j in rng:
            for k in rng:
                R[i][j] += A[i][k] * B[k][j]
    return R

print mult(M, M_inv)

[[1, 0, 0], [0, 1, 0], [0, 0, 1]]