Implement decomposition based solve functions for solving Ax = b

This is a huge help in debugging. So I decided to keep it. We can remove this later.
How do we make the following work when A is a matrix.

```C++
std::cout << A << std::endl
```
I have to use `A.__str__()` to print the matrix into output stream.