-// superfluous helper function
-ex matrix::ffe_get(unsigned r, unsigned c) const
-{
- return operator()(r-1,c-1);
-}
-
-/** Solve a set of equations for an m x n matrix by fraction-free Gaussian
- * elimination. Based on algorithm 9.1 from 'Algorithms for Computer Algebra'
- * by Keith O. Geddes et al.
- *
- * @param vars n x p matrix
- * @param rhs m x p matrix
- * @exception logic_error (incompatible matrices)
- * @exception runtime_error (singular matrix) */
-matrix matrix::fraction_free_elim(const matrix & vars,
- const matrix & rhs) const
-{
- // FIXME: implement a Sasaki-Murao scheme which avoids division at all!
- if ((row != rhs.row) || (col != vars.row) || (rhs.col != vars.col))
- throw (std::logic_error("matrix::fraction_free_elim(): incompatible matrices"));
-
- matrix a(*this); // make a copy of the matrix
- matrix b(rhs); // make a copy of the rhs vector
-
- // given an m x n matrix a, reduce it to upper echelon form
- unsigned m = a.row;
- unsigned n = a.col;
- int sign = 1;
- ex divisor = 1;
- unsigned r = 1;
-
- // eliminate below row r, with pivot in column k
- for (unsigned k=1; (k<=n)&&(r<=m); ++k) {
- // find a nonzero pivot
- unsigned p;
- for (p=r; (p<=m)&&(a.ffe_get(p,k).is_equal(_ex0())); ++p) {}
- // pivot is in row p
- if (p<=m) {
- if (p!=r) {
- // switch rows p and r
- for (unsigned j=k; j<=n; ++j)
- a.ffe_swap(p,j,r,j);
- b.ffe_swap(p,1,r,1);
- // keep track of sign changes due to row exchange
- sign = -sign;
- }
- for (unsigned i=r+1; i<=m; ++i) {
- for (unsigned j=k+1; j<=n; ++j) {
- a.ffe_set(i,j,(a.ffe_get(r,k)*a.ffe_get(i,j)
- -a.ffe_get(r,j)*a.ffe_get(i,k))/divisor);
- a.ffe_set(i,j,a.ffe_get(i,j).normal() /*.normal() */ );
- }
- b.ffe_set(i,1,(a.ffe_get(r,k)*b.ffe_get(i,1)
- -b.ffe_get(r,1)*a.ffe_get(i,k))/divisor);
- b.ffe_set(i,1,b.ffe_get(i,1).normal() /*.normal() */ );
- a.ffe_set(i,k,0);
- }
- divisor = a.ffe_get(r,k);
- r++;
- }
- }
- // optionally compute the determinant for square or augmented matrices
- // if (r==m+1) { det = sign*divisor; } else { det = 0; }
-
- /*
- for (unsigned r=1; r<=m; ++r) {
- for (unsigned c=1; c<=n; ++c) {
- cout << a.ffe_get(r,c) << "\t";
- }
- cout << " | " << b.ffe_get(r,1) << endl;
- }
- */
-
-#ifdef DO_GINAC_ASSERT
- // test if we really have an upper echelon matrix
- int zero_in_last_row = -1;
- for (unsigned r=1; r<=m; ++r) {
- int zero_in_this_row=0;
- for (unsigned c=1; c<=n; ++c) {
- if (a.ffe_get(r,c).is_equal(_ex0()))
- zero_in_this_row++;
- else
- break;
- }
- GINAC_ASSERT((zero_in_this_row>zero_in_last_row)||(zero_in_this_row=n));
- zero_in_last_row = zero_in_this_row;
- }
-#endif // def DO_GINAC_ASSERT
-
- /*
- cout << "after" << endl;
- cout << "a=" << a << endl;
- cout << "b=" << b << endl;
- */
-
- // assemble solution
- matrix sol(n,1);
- unsigned last_assigned_sol = n+1;
- for (unsigned r=m; r>0; --r) {
- unsigned first_non_zero = 1;
- while ((first_non_zero<=n)&&(a.ffe_get(r,first_non_zero).is_zero()))
- first_non_zero++;
- if (first_non_zero>n) {
- // row consists only of zeroes, corresponding rhs must be 0 as well
- if (!b.ffe_get(r,1).is_zero()) {
- throw (std::runtime_error("matrix::fraction_free_elim(): singular matrix"));
- }
- } else {
- // assign solutions for vars between first_non_zero+1 and
- // last_assigned_sol-1: free parameters
- for (unsigned c=first_non_zero+1; c<=last_assigned_sol-1; ++c) {
- sol.ffe_set(c,1,vars.ffe_get(c,1));
- }
- ex e = b.ffe_get(r,1);
- for (unsigned c=first_non_zero+1; c<=n; ++c) {
- e=e-a.ffe_get(r,c)*sol.ffe_get(c,1);
- }
- sol.ffe_set(first_non_zero,1,
- (e/a.ffe_get(r,first_non_zero)).normal());
- last_assigned_sol = first_non_zero;
- }
- }
- // assign solutions for vars between 1 and
- // last_assigned_sol-1: free parameters
- for (unsigned c=1; c<=last_assigned_sol-1; ++c)
- sol.ffe_set(c,1,vars.ffe_get(c,1));
-
-#ifdef DO_GINAC_ASSERT
- // test solution with echelon matrix
- for (unsigned r=1; r<=m; ++r) {
- ex e = 0;
- for (unsigned c=1; c<=n; ++c)
- e = e+a.ffe_get(r,c)*sol.ffe_get(c,1);
- if (!(e-b.ffe_get(r,1)).normal().is_zero()) {
- cout << "e=" << e;
- cout << "b.ffe_get(" << r<<",1)=" << b.ffe_get(r,1) << endl;
- cout << "diff=" << (e-b.ffe_get(r,1)).normal() << endl;
- }
- GINAC_ASSERT((e-b.ffe_get(r,1)).normal().is_zero());
- }
-
- // test solution with original matrix
- for (unsigned r=1; r<=m; ++r) {
- ex e = 0;
- for (unsigned c=1; c<=n; ++c)
- e = e+ffe_get(r,c)*sol.ffe_get(c,1);
- try {
- if (!(e-rhs.ffe_get(r,1)).normal().is_zero()) {
- cout << "e=" << e << endl;
- e.printtree(cout);
- ex en = e.normal();
- cout << "e.normal()=" << en << endl;
- en.printtree(cout);
- cout << "rhs.ffe_get(" << r<<",1)=" << rhs.ffe_get(r,1) << endl;
- cout << "diff=" << (e-rhs.ffe_get(r,1)).normal() << endl;
- }
- } catch (...) {
- ex xxx = e - rhs.ffe_get(r,1);
- cerr << "xxx=" << xxx << endl << endl;
- }
- GINAC_ASSERT((e-rhs.ffe_get(r,1)).normal().is_zero());
- }
-#endif // def DO_GINAC_ASSERT
-
- return sol;
-}