I feel this matrix class is an improvement from the previous one, but it still lacks the robustness I need. Moreover, there are still some room for performance improvements and new functions. Here are below the ideas I have to improve this matrix class :
1- Reoptimise code by using the convention that "zero" is not a stored in the vectors and thus, replace for by foreach in most of the routines.
2- Introduce a general Gauss-Seidel algorithm that is more numerically stable than "classic" matrix inversion. This will probably yield new functions such diag (the diagonal of a matrix), upper_tri (upper triangle matrix), and lower_tri (lower triangle matrix) for matrix decomposition.
3- Introduce a method to find eigenvalues (and eigenvectors).
4- Some basic functions such as the trace of a matrix and whatever else that cross my mind which does not take 2000 lines of code.
More generally, I intend to work on :
5- Matrix representation of quadratic approximation of a system of equations and Quadratic regulators.
6- Automatised numerical estimation of the previous point based on a set of "code like" set of equations (if necessary).