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

[Vitaly Magerya]

Hi, folks. I've run into a number of problems with matrix solution
routines in GiNaC, and I'd like to sort them out one by one.

The basic overview of the problems is this:

1) matrix::gauss_elimination checks for zero entries via a
   simple .is_zero() call without a .normal() call before, so it
   misses zeros in non-normal form. For example it will not detect
   that this expression:

       1/(x - 1) + 1/(x + 1) - 2*x/(x*x - 1)

   ... is zero. This leads to "division by zero" errors.

   Test case:

       #include <ginac/ginac.h>
       using namespace GiNaC;
       int main() {
           symbol x("x");
           matrix mx {
               {1,0,0},
               {0,1/(x+1)-1/(x-1)+2/(x*x-1),1},
               {0,1,2},
           };
           matrix zero(mx.rows(), 1);
           matrix tmp(mx.cols(), 1);
           for (unsigned i = 0; i < tmp.nops(); i++) {
               tmp.let_op(i) = symbol();
           }
           mx.solve(tmp, zero, solve_algo::gauss);
           return 0;
       }

2) matrix matrix::solve does the same erorr, so if gauss_elimination
   returned a matrix with un-normal zeros, it too may error out
   with "division by zero".

   Test case:

       #include <ginac/ginac.h>
       using namespace GiNaC;
       int main() {
           symbol x("x");
           matrix mx {
               {1,0,0},
               {0,1/(x+1)-1/(x-1)+2/(x*x-1),1},
               {0,0,0},
           };
           matrix zero(mx.rows(), 1);
           matrix tmp(mx.cols(), 1);
           for (unsigned i = 0; i < tmp.nops(); i++) {
               tmp.let_op(i) = symbol();
           }
           mx.solve(tmp, zero, solve_algo::gauss);
           return 0;
       }

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.

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(): it doesn't provide an 'algo'
argument, and always defaults to solve_algo::automatic (which
means fraction_free_elimination). So for the start I propose to
add that argument (see the attached patch). The reservation I have
with this patch is other functions (solve() and determinant())
have their own enumerations for the 'algo' argument, while
inverse() will be reusing the one from solve(). Not the cleanest
thing in the world from API usage point of view. What do you
think of it?
-------------- next part --------------
diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index 915d4543..14c27d50 100644
--- a/ginac/matrix.cpp
+++ b/ginac/matrix.cpp
@@ -942,7 +942,7 @@ ex matrix::charpoly(const ex & lambda) const
  *  @return    the inverted matrix
  *  @exception logic_error (matrix not square)
  *  @exception runtime_error (singular matrix) */
-matrix matrix::inverse() const
+matrix matrix::inverse(unsigned algo) const
 {
 	if (row != col)
 		throw (std::logic_error("matrix::inverse(): matrix not square"));
@@ -965,7 +965,7 @@ matrix matrix::inverse() const
 	
 	matrix sol(row,col);
 	try {
-		sol = this->solve(vars,identity);
+		sol = this->solve(vars,identity,algo);
 	} catch (const std::runtime_error & e) {
 	    if (e.what()==std::string("matrix::solve(): inconsistent linear system"))
 			throw (std::runtime_error("matrix::inverse(): singular matrix"));
diff --git a/ginac/matrix.h b/ginac/matrix.h
index 1bb549b0..ce26d163 100644
--- a/ginac/matrix.h
+++ b/ginac/matrix.h
@@ -148,7 +148,7 @@ public:
 	ex determinant(unsigned algo = determinant_algo::automatic) const;
 	ex trace() const;
 	ex charpoly(const ex & lambda) const;
-	matrix inverse() const;
+	matrix inverse(unsigned algo = solve_algo::automatic) const;
 	matrix solve(const matrix & vars, const matrix & rhs,
 	             unsigned algo = solve_algo::automatic) const;
 	unsigned rank() const;