Implement decomposition based solve functions for solving Ax = b

You should keep this, we need good Matrix printing. In order to hook it up into the stream, have a look at the bottom of `dict.h`:
```
std::ostream& operator<<(std::ostream& out, const CSymPy::umap_basic_int& d);
```
You need to create a similar function for `const DenseMatrix &`, then it should work.