Implement decomposition based solve functions for solving Ax = b

I changed the API of this function to return one Matrix which contain both `L` and `U`. Here `L` is the lower triangular part of `A` and `U` is the upper triangular part. This make sense as we don't need to create two matrices when it comes to solving the system `Ax = b`.