- * @exception runtime_error (singular matrix) */
-matrix matrix::fraction_free_elim(const matrix & vars,
- const matrix & rhs) const
-{
- // FIXME: use implementation of matrix::fraction_free_elimination instead!
- 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 = 0;
-
- // eliminate below row r, with pivot in column k
- for (unsigned k=0; (k<n)&&(r<m); ++k) {
- // find a nonzero pivot
- unsigned p;
- for (p=r; (p<m)&&(a.m[p*a.cols()+k].is_zero()); ++p) {}
- // pivot is in row p
- if (p<m) {
- if (p!=r) {
- // swap rows p and r
- for (unsigned j=k; j<n; ++j)
- a.m[p*a.cols()+j].swap(a.m[r*a.cols()+j]);
- b.m[p*b.cols()].swap(b.m[r*b.cols()]);
- // keep track of sign changes due to row exchange
- sign *= -1;
- }
- for (unsigned i=r+1; i<m; ++i) {
- for (unsigned j=k+1; j<n; ++j) {
- a.set(i,j,(a.m[r*a.cols()+k]*a.m[i*a.cols()+j]
- -a.m[r*a.cols()+j]*a.m[i*a.cols()+k])/divisor);
- a.set(i,j,a.m[i*a.cols()+j].normal());
- }
- b.set(i,0,(a.m[r*a.cols()+k]*b.m[i*b.cols()]
- -b.m[r*b.cols()]*a.m[i*a.cols()+k])/divisor);
- b.set(i,0,b.m[i*b.cols()].normal());
- a.set(i,k,_ex0());
- }
- divisor = a.m[r*a.cols()+k];
- ++r;
- }
- }
-
-#ifdef DO_GINAC_ASSERT
- // test if we really have an upper echelon matrix
- int zero_in_last_row = -1;
- for (unsigned r=0; r<m; ++r) {
- int zero_in_this_row=0;
- for (unsigned c=0; c<n; ++c) {
- if (a.m[r*a.cols()+c].is_zero())
- ++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
-
- // assemble solution
- matrix sol(n,1);
- unsigned last_assigned_sol = n+1;
- for (int r=m-1; r>=0; --r) {
- unsigned first_non_zero = 1;
- while ((first_non_zero<=n)&&(a(r,first_non_zero-1).is_zero()))
- first_non_zero++;
- if (first_non_zero>n) {
- // row consists only of zeroes, corresponding rhs must be 0 as well
- if (!b.m[r*b.cols()].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; c<last_assigned_sol-1; ++c)
- sol.set(c,0,vars.m[c*vars.cols()]);
- ex e = b.m[r*b.cols()];
- for (unsigned c=first_non_zero; c<n; ++c)
- e -= a.m[r*a.cols()+c]*sol.m[c*sol.cols()];
- sol.set(first_non_zero-1,0,
- (e/(a.m[r*a.cols()+(first_non_zero-1)])).normal());
- last_assigned_sol = first_non_zero;
- }
- }
- // assign solutions for vars between 1 and
- // last_assigned_sol-1: free parameters
- for (unsigned c=0; c<last_assigned_sol-1; ++c)
- sol.set(c,0,vars.m[c*vars.cols()]);
-
-#ifdef DO_GINAC_ASSERT
- // test solution with echelon matrix
- for (unsigned r=0; r<m; ++r) {
- ex e = 0;
- for (unsigned c=0; c<n; ++c)
- e += a(r,c)*sol(c,0);
- if (!(e-b(r,0)).normal().is_zero()) {
- cout << "e=" << e;
- cout << "b(" << r <<",0)=" << b(r,0) << endl;
- cout << "diff=" << (e-b(r,0)).normal() << endl;
- }
- GINAC_ASSERT((e-b(r,0)).normal().is_zero());
- }
-
- // test solution with original matrix
- for (unsigned r=0; r<m; ++r) {
- ex e = 0;
- for (unsigned c=0; c<n; ++c)
- e += this->m[r*cols()+c]*sol(c,0);
- try {
- if (!(e-rhs(r,0)).normal().is_zero()) {
- cout << "e==" << e << endl;
- e.printtree(cout);
- ex en = e.normal();
- cout << "e.normal()=" << en << endl;
- en.printtree(cout);
- cout << "rhs(" << r <<",0)=" << rhs(r,0) << endl;
- cout << "diff=" << (e-rhs(r,0)).normal() << endl;
- }
- } catch (...) {
- ex xxx = e - rhs(r,0);
- cerr << "xxx=" << xxx << endl << endl;
- }
- GINAC_ASSERT((e-rhs(r,0)).normal().is_zero());
- }
-#endif // def DO_GINAC_ASSERT
-
- return sol;
-}
-
-/** Solve a set of equations for an m x n matrix.
- *
- * @param vars n x p matrix
- * @param rhs m x p matrix
- * @exception logic_error (incompatible matrices)
- * @exception runtime_error (singular matrix) */