Implement decomposition based solve functions for solving Ax = b

This should be
```C++
void fraction_free_LU_solve(RCP<const DenseMatrix> &A, RCP<const DenseMatrix> &b, Ptr<RCP<DenseMatrix>> &x)
```
rather than the current definition. But we can't use `RCP` with matrix classes as we don't derive it from `Basic`.

I guess there is no harm in deriving from `Basic`. We can have some current APIs as they are now and can better address these situations. Also, say `CSymPy` is going to implement a `set` class in future which models a mathematical set, then we will need to have matrices in it as well.