On 01/18/2018 03:10 PM, I wrote:
[...] Now, there are multiple ways of dealing with these problems. The practical workaround for the moment is to: a) use solve_algo::gauss for all sparse-ish matrices b) normal()ize the whole matrix before passing it into solve()
There are places where this workaround is not available. One of such places is matrix::inverse()
The second place where a better algorithm is needed is matrix::rank(). Currently it blindly fraction_free_elimination(), which should never be used with sparse matrices. BTW, I noticed that matrix::determinant() has some logic to determine matrix sparseness, and it would actually switch away from fraction_free_elimination() for sparse matrices (it will use minor expansion method instead), where sparseness is determined by comparing the number of zero cells to the number of all cells divided by 5. I think it would be beneficial if the same sort of logic would be used in matrix::solve() too. (See the attached patch for my proposal). This would additionally alleviate the need for an 'algo' argument in matrix::inverse(), since solve_algo::automatic would now perform the sensible thing for sparse matrices (which is Gauss elimination). This of course doesn't solve the problem with matrix::rank(). I think there are two good ways to proceed: 1) Apply the attached patch, and make matrix::rank() use similar logic to determine which *_elimination() to call. 2) Apply the attached patch, and add 'algo' arguments to both matrix::inverse() and matrix::rank() (probably using separate enumerations for each).