Hi, folks. Since we've recently added solve_algo::markowitz, I thought it would be a good idea to update the automatic elimination algorithm selection logic we have for matrix::solve() (currently living in matrix::echelon_form()). To this end I've constructed a benchmark that generates random-ish matrices and measures how fast matrix::solve() can solve them. The matrix generation works by starting with a known simple solution and then permuting and mixing rows and columns randomly. You can find the code for this procedure at [1]. With that code I've generated a bunch of data points and run a quick analysis on it. Overall, the fastest (on average) algorithm seems to be: For numeric matrices: If ncells > 200 && density < 0.5: Markowitz Otherwise: Gauss For symbolic matrices: If density > 0.6 && ncells <= 12: Divfree If density > 0.6 && ncells < 120: Bareiss Otherwise: Markowitz You can find the details (and some pretty charts) at [2]. The attached patch implements this recommendation. Note that matrix::determinant() also has a similar automatic algorithm selection logic, but I didn't touch it for now. [1] https://github.com/magv/ginac-matrix-solve-bench/blob/master/mx-solve-bench.... [2] https://github.com/magv/ginac-matrix-solve-bench/blob/master/analysis.ipynb