import numpy as np

coefficients = np.array(
    (
        (3.0, 1.0, 2.0),
        (2.0, 6.0, 1.0),
        (1.0, -1.0, 4.0),
    )
)
rhs = np.array((11.0, 1.0, 15.0))

solution = np.linalg.solve(coefficients, rhs)
reconstructed = coefficients @ solution
residual = reconstructed - rhs

print("solution:", np.round(solution, 6))
print("reconstructed rhs:", np.round(reconstructed, 6))
print("residual norm:", round(float(np.linalg.norm(residual)), 12))
print("matches rhs:", np.allclose(reconstructed, rhs))