Implement decomposition based solve functions for solving Ax = b

That's unrelated to this PR, but I think it's a better design to just call some `eq()` method from the matrix class and just call that method from here, i.e. just quickly hook it up into the C++ `==` operation.