On 18.01.2018 15:10, Vitaly Magerya wrote:
3) The default echelon form algorithm, Bareiss elimination, (matrix::fraction_free_elimination) is unbearably slow for sparse matrices, and should never be used with them. This is because of explosive common factor growth during its evaluation.
For an illustration, using matrix::fraction_free_elimination on a matrix like this:
{{a11,a12,a13,a14}, { 0,a22,a23,a24}, { 0, 0,a33,a34}, { 0, 0, 0,a44}}
... will give you a result like this:
{{a11, a12, a13, a14}, { 0,a11*a22, a11*a23, a11*a24}, { 0, 0,a11*a22*a33, a11*a22*a34}, { 0, 0, 0,a11*a22*a33*a44}}
Now imagine a similar 100x100 matrix with slightly more complex diagonal elements.
To consider an echelon input matrix is indeed an instructive example. In the first elimination step, Bareiss fraction free elimination produces a 3x3 submatrix where all elements have the common factor a11: {{a11, a12, a13, a14}, { 0,a11*a22,a11*a23,a11*a24}, { 0, 0,a11*a33,a11*a34}, { 0, 0, 0,a11*a44}} This is because the elements in the first column below the pivot element a11 are all zero. It should be possible to use this observation and simply skip this fraction free elimination step! This would avoid the term growth if the elements below the pivot are all zero. Of course, we mustn't divide by a11 in the next step (because it won't divide). And this won't work for determinant computation (because the determinant ends up in the lower right element with the Bareiss elimination). I'm not against tuning the heuristic choice which algorithm to use. But with this modified fraction free elimination, the input matrix would be equal to the output matrix: No term growth for matrices in echelon form and no GCD computation ever. Cheers -richy.