[GiNaC-list] matrix::solve()-related problems

[Richard B. Kreckel]

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.