- if ((row != rhs.row) || (col != vars.row) || (rhs.col != vars.col)) {
- throw (std::logic_error("matrix::solve(): 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
- int m=a.row;
- int n=a.col;
- int sign=1;
- ex divisor=1;
- int r=1;
-
- // eliminate below row r, with pivot in column k
- for (int k=1; (k<=n)&&(r<=m); ++k) {
- // find a nonzero pivot
- int p;
- for (p=r; (p<=m)&&(a.ffe_get(p,k).is_equal(exZERO())); ++p) {}
- // pivot is in row p
- if (p<=m) {
- if (p!=r) {
- // switch rows p and r
- for (int 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 (int i=r+1; i<=m; ++i) {
- for (int 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 (int r=1; r<=m; ++r) {
- for (int 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 (int r=1; r<=m; ++r) {
- int zero_in_this_row=0;
- for (int c=1; c<=n; ++c) {
- if (a.ffe_get(r,c).is_equal(exZERO())) {
- 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);
- int last_assigned_sol=n+1;
- for (int r=m; r>0; --r) {
- int 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 (int 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 (int 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 (int c=1; c<=last_assigned_sol-1; ++c) {
- sol.ffe_set(c,1,vars.ffe_get(c,1));
- }
-
- /*
- for (int c=1; c<=n; ++c) {
- cout << vars.ffe_get(c,1) << "->" << sol.ffe_get(c,1) << endl;
- }
- */
-
-#ifdef DO_GINAC_ASSERT
- // test solution with echelon matrix
- for (int r=1; r<=m; ++r) {
- ex e=0;
- for (int 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 (int r=1; r<=m; ++r) {
- ex e=0;
- for (int 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;
-}
-
-/** Solve simultaneous set of equations. */
-matrix matrix::solve(matrix const & v) const
+ ensure_if_modifiable();
+ const unsigned m = this->rows();
+ const unsigned n = this->cols();
+ GINAC_ASSERT(!det || n==m);
+ int sign = 1;
+
+ unsigned r0 = 0;
+ for (unsigned r1=0; (r1<n-1)&&(r0<m-1); ++r1) {
+ int indx = pivot(r0, r1, true);
+ if (indx == -1) {
+ sign = 0;
+ if (det)
+ return 0; // leaves *this in a messy state
+ }
+ if (indx>=0) {
+ if (indx > 0)
+ sign = -sign;
+ for (unsigned r2=r0+1; r2<m; ++r2) {
+ if (!this->m[r2*n+r1].is_zero()) {
+ // yes, there is something to do in this row
+ ex piv = this->m[r2*n+r1] / this->m[r0*n+r1];
+ for (unsigned c=r1+1; c<n; ++c) {
+ this->m[r2*n+c] -= piv * this->m[r0*n+c];
+ if (!this->m[r2*n+c].info(info_flags::numeric))
+ this->m[r2*n+c] = this->m[r2*n+c].normal();
+ }
+ }
+ // fill up left hand side with zeros
+ for (unsigned c=0; c<=r1; ++c)
+ this->m[r2*n+c] = _ex0();
+ }
+ if (det) {
+ // save space by deleting no longer needed elements
+ for (unsigned c=r0+1; c<n; ++c)
+ this->m[r0*n+c] = _ex0();
+ }
+ ++r0;
+ }
+ }
+
+ return sign;
+}
+
+
+/** Perform the steps of division free elimination to bring the m x n matrix
+ * into an upper echelon form.
+ *
+ * @param det may be set to true to save a lot of space if one is only
+ * interested in the diagonal elements (i.e. for calculating determinants).
+ * The others are set to zero in this case.
+ * @return sign is 1 if an even number of rows was swapped, -1 if an odd
+ * number of rows was swapped and 0 if the matrix is singular. */
+int matrix::division_free_elimination(const bool det)