- Added method matrix::rank().

- Fixes to matrix elimination routines, needed for rank computations.

diff --git a/check/exam_matrices.cpp b/check/exam_matrices.cpp
index 67c52be85b841c36f15f86f27e1a48b7b50df61f..449b6d38a15b6f47c280d823156d507b77c4b434 100644 (file)
--- a/check/exam_matrices.cpp
+++ b/check/exam_matrices.cpp
@@ -241,6 +241,46 @@ static unsigned matrix_evalm()
        return result;
 }
 
+static unsigned matrix_rank()
+{
+       unsigned result = 0;
+       symbol x("x"), y("y");
+       matrix m(3,3);
+
+       // the zero matrix always has rank 0
+       if (m.rank() != 0) {
+               clog << "The rank of " << m << " was not computed correctly." << endl;
+               ++result;
+       }
+
+       // a trivial rank one example
+       m = 1, 0, 0,
+           2, 0, 0,
+           3, 0, 0;
+       if (m.rank() != 1) {
+               clog << "The rank of " << m << " was not computed correctly." << endl;
+               ++result;
+       }
+
+       // an example from Maple's help with rank two
+       m = x,  1,  0,
+           0,  0,  1,
+          x*y, y,  1;
+       if (m.rank() != 2) {
+               clog << "The rank of " << m << " was not computed correctly." << endl;
+               ++result;
+       }
+
+       // the 3x3 unit matrix has rank 3
+       m = ex_to<matrix>(unit_matrix(3,3));
+       if (m.rank() != 3) {
+               clog << "The rank of " << m << " was not computed correctly." << endl;
+               ++result;
+       }
+
+       return result;  
+}
+
 static unsigned matrix_misc()
 {
        unsigned result = 0;
@@ -305,6 +345,7 @@ unsigned exam_matrices()
        result += matrix_invert3();  cout << '.' << flush;
        result += matrix_solve2();  cout << '.' << flush;
        result += matrix_evalm();  cout << "." << flush;
+       result += matrix_rank();  cout << "." << flush;
        result += matrix_misc();  cout << '.' << flush;
        
        if (!result) {

diff --git a/doc/tutorial/ginac.texi b/doc/tutorial/ginac.texi
index 3f98183162437c3a8d7c4bcf03c9cc977c52e2eb..b58d631b868143183d24b5777441e7d17e93ae61 100644 (file)
--- a/doc/tutorial/ginac.texi
+++ b/doc/tutorial/ginac.texi
@@ -2013,15 +2013,17 @@ more information about using matrices with indices, and about indices in
 general.
 
 The @code{matrix} class provides a couple of additional methods for
-computing determinants, traces, and characteristic polynomials:
+computing determinants, traces, characteristic polynomials and ranks:
 
 @cindex @code{determinant()}
 @cindex @code{trace()}
 @cindex @code{charpoly()}
+@cindex @code{rank()}
 @example
 ex matrix::determinant(unsigned algo=determinant_algo::automatic) const;
 ex matrix::trace() const;
 ex matrix::charpoly(const ex & lambda) const;
+unsigned matrix::rank() const;
 @end example
 
 The @samp{algo} argument of @code{determinant()} allows to select

diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index 2c6e5bb4406f9559f36462c9df3038543e8be1e3..e24f644741d6be1745c610020d314b5c338cf7ac 100644 (file)
--- a/ginac/matrix.cpp
+++ b/ginac/matrix.cpp
@@ -764,7 +764,7 @@ ex matrix::determinant(unsigned algo) const
                else
                        return m[0].expand();
        }
-       
+
        // Compute the determinant
        switch(algo) {
                case determinant_algo::gauss: {
@@ -860,7 +860,7 @@ ex matrix::trace() const
                tr += m[r*col+r];
        
        if (tr.info(info_flags::rational_function) &&
-               !tr.info(info_flags::crational_polynomial))
+          !tr.info(info_flags::crational_polynomial))
                return tr.normal();
        else
                return tr.expand();
@@ -973,8 +973,8 @@ matrix matrix::inverse() const
  *  @exception runtime_error (inconsistent linear system)
  *  @see       solve_algo */
 matrix matrix::solve(const matrix & vars,
-                                        const matrix & rhs,
-                                        unsigned algo) const
+                     const matrix & rhs,
+                     unsigned algo) const
 {
        const unsigned m = this->rows();
        const unsigned n = this->cols();
@@ -1067,6 +1067,29 @@ matrix matrix::solve(const matrix & vars,
 }
 
 
+/** Compute the rank of this matrix. */
+unsigned matrix::rank() const
+{
+       // Method:
+       // Transform this matrix into upper echelon form and then count the
+       // number of non-zero rows.
+
+       GINAC_ASSERT(row*col==m.capacity());
+
+       // Actually, any elimination scheme will do since we are only
+       // interested in the echelon matrix' zeros.
+       matrix to_eliminate = *this;
+       to_eliminate.fraction_free_elimination();
+
+       unsigned r = row*col;  // index of last non-zero element
+       while (r--) {
+               if (!to_eliminate.m[r].is_zero())
+                       return 1+r/col;
+       }
+       return 0;
+}
+
+
 // protected
 
 /** Recursive determinant for small matrices having at least one symbolic
@@ -1208,8 +1231,8 @@ int matrix::gauss_elimination(const bool det)
        int sign = 1;
        
        unsigned r0 = 0;
-       for (unsigned r1=0; (r1<n-1)&&(r0<m-1); ++r1) {
-               int indx = pivot(r0, r1, true);
+       for (unsigned c0=0; c0<n && r0<m-1; ++c0) {
+               int indx = pivot(r0, c0, true);
                if (indx == -1) {
                        sign = 0;
                        if (det)
@@ -1219,17 +1242,17 @@ int matrix::gauss_elimination(const bool det)
                        if (indx > 0)
                                sign = -sign;
                        for (unsigned r2=r0+1; r2<m; ++r2) {
-                               if (!this->m[r2*n+r1].is_zero()) {
+                               if (!this->m[r2*n+c0].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) {
+                                       ex piv = this->m[r2*n+c0] / this->m[r0*n+c0];
+                                       for (unsigned c=c0+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)
+                               for (unsigned c=r0; c<=c0; ++c)
                                        this->m[r2*n+c] = _ex0;
                        }
                        if (det) {
@@ -1240,7 +1263,12 @@ int matrix::gauss_elimination(const bool det)
                        ++r0;
                }
        }
-       
+       // clear remaining rows
+       for (unsigned r=r0+1; r<m; ++r) {
+               for (unsigned c=0; c<n; ++c)
+                       this->m[r*n+c] = _ex0;
+       }
+
        return sign;
 }
 
@@ -1262,8 +1290,8 @@ int matrix::division_free_elimination(const bool det)
        int sign = 1;
        
        unsigned r0 = 0;
-       for (unsigned r1=0; (r1<n-1)&&(r0<m-1); ++r1) {
-               int indx = pivot(r0, r1, true);
+       for (unsigned c0=0; c0<n && r0<m-1; ++c0) {
+               int indx = pivot(r0, c0, true);
                if (indx==-1) {
                        sign = 0;
                        if (det)
@@ -1273,10 +1301,10 @@ int matrix::division_free_elimination(const bool det)
                        if (indx>0)
                                sign = -sign;
                        for (unsigned r2=r0+1; r2<m; ++r2) {
-                               for (unsigned c=r1+1; c<n; ++c)
-                                       this->m[r2*n+c] = (this->m[r0*n+r1]*this->m[r2*n+c] - this->m[r2*n+r1]*this->m[r0*n+c]).expand();
+                               for (unsigned c=c0+1; c<n; ++c)
+                                       this->m[r2*n+c] = (this->m[r0*n+c0]*this->m[r2*n+c] - this->m[r2*n+c0]*this->m[r0*n+c]).expand();
                                // fill up left hand side with zeros
-                               for (unsigned c=0; c<=r1; ++c)
+                               for (unsigned c=r0; c<=c0; ++c)
                                        this->m[r2*n+c] = _ex0;
                        }
                        if (det) {
@@ -1287,7 +1315,12 @@ int matrix::division_free_elimination(const bool det)
                        ++r0;
                }
        }
-       
+       // clear remaining rows
+       for (unsigned r=r0+1; r<m; ++r) {
+               for (unsigned c=0; c<n; ++c)
+                       this->m[r*n+c] = _ex0;
+       }
+
        return sign;
 }
 
@@ -1360,8 +1393,8 @@ int matrix::fraction_free_elimination(const bool det)
        }
        
        unsigned r0 = 0;
-       for (unsigned r1=0; (r1<n-1)&&(r0<m-1); ++r1) {
-               int indx = tmp_n.pivot(r0, r1, true);
+       for (unsigned c0=0; c0<n && r0<m-1; ++c0) {
+               int indx = tmp_n.pivot(r0, c0, true);
                if (indx==-1) {
                        sign = 0;
                        if (det)
@@ -1371,17 +1404,17 @@ int matrix::fraction_free_elimination(const bool det)
                        if (indx>0) {
                                sign = -sign;
                                // tmp_n's rows r0 and indx were swapped, do the same in tmp_d:
-                               for (unsigned c=r1; c<n; ++c)
+                               for (unsigned c=c0; c<n; ++c)
                                        tmp_d.m[n*indx+c].swap(tmp_d.m[n*r0+c]);
                        }
                        for (unsigned r2=r0+1; r2<m; ++r2) {
-                               for (unsigned c=r1+1; c<n; ++c) {
-                                       dividend_n = (tmp_n.m[r0*n+r1]*tmp_n.m[r2*n+c]*
-                                                     tmp_d.m[r2*n+r1]*tmp_d.m[r0*n+c]
-                                                    -tmp_n.m[r2*n+r1]*tmp_n.m[r0*n+c]*
-                                                     tmp_d.m[r0*n+r1]*tmp_d.m[r2*n+c]).expand();
-                                       dividend_d = (tmp_d.m[r2*n+r1]*tmp_d.m[r0*n+c]*
-                                                     tmp_d.m[r0*n+r1]*tmp_d.m[r2*n+c]).expand();
+                               for (unsigned c=c0+1; c<n; ++c) {
+                                       dividend_n = (tmp_n.m[r0*n+c0]*tmp_n.m[r2*n+c]*
+                                                     tmp_d.m[r2*n+c0]*tmp_d.m[r0*n+c]
+                                                    -tmp_n.m[r2*n+c0]*tmp_n.m[r0*n+c]*
+                                                     tmp_d.m[r0*n+c0]*tmp_d.m[r2*n+c]).expand();
+                                       dividend_d = (tmp_d.m[r2*n+c0]*tmp_d.m[r0*n+c]*
+                                                     tmp_d.m[r0*n+c0]*tmp_d.m[r2*n+c]).expand();
                                        bool check = divide(dividend_n, divisor_n,
                                                            tmp_n.m[r2*n+c], true);
                                        check &= divide(dividend_d, divisor_d,
@@ -1389,13 +1422,13 @@ int matrix::fraction_free_elimination(const bool det)
                                        GINAC_ASSERT(check);
                                }
                                // fill up left hand side with zeros
-                               for (unsigned c=0; c<=r1; ++c)
+                               for (unsigned c=r0; c<=c0; ++c)
                                        tmp_n.m[r2*n+c] = _ex0;
                        }
-                       if ((r1<n-1)&&(r0<m-1)) {
+                       if (c0<n && r0<m-1) {
                                // compute next iteration's divisor
-                               divisor_n = tmp_n.m[r0*n+r1].expand();
-                               divisor_d = tmp_d.m[r0*n+r1].expand();
+                               divisor_n = tmp_n.m[r0*n+c0].expand();
+                               divisor_d = tmp_d.m[r0*n+c0].expand();
                                if (det) {
                                        // save space by deleting no longer needed elements
                                        for (unsigned c=0; c<n; ++c) {
@@ -1407,6 +1440,12 @@ int matrix::fraction_free_elimination(const bool det)
                        ++r0;
                }
        }
+       // clear remaining rows
+       for (unsigned r=r0+1; r<m; ++r) {
+               for (unsigned c=0; c<n; ++c)
+                       tmp_n.m[r*n+c] = _ex0;
+       }
+
        // repopulate *this matrix:
        exvector::iterator it = this->m.begin(), itend = this->m.end();
        tmp_n_it = tmp_n.m.begin();

diff --git a/ginac/matrix.h b/ginac/matrix.h
index ee30601310a2a7d8646bb51d9ef3410a11bb2e57..29af2a0747d903198fa1cfa78347a5e2ab0097ba 100644 (file)
--- a/ginac/matrix.h
+++ b/ginac/matrix.h
@@ -148,6 +148,7 @@ public:
        matrix inverse() const;
        matrix solve(const matrix & vars, const matrix & rhs,
                     unsigned algo = solve_algo::automatic) const;
+       unsigned rank() const;
 protected:
        ex determinant_minor() const;
        int gauss_elimination(const bool det = false);
@@ -203,6 +204,9 @@ inline ex charpoly(const matrix & m, const ex & lambda)
 inline matrix inverse(const matrix & m)
 { return m.inverse(); }
 
+inline unsigned rank(const matrix & m)
+{ return m.rank(); }
+
 // utility functions
 
 /** Specialization of is_exactly_a<matrix>(obj) for matrix objects. */

diff --git a/ginsh/ginsh.1.in b/ginsh/ginsh.1.in
index 01ecbb59a46b71188ee96e93dfe4270601c1a89d..1662ce2922f71b8dc813c4e6282acc78487d10c7 100644 (file)
--- a/ginsh/ginsh.1.in
+++ b/ginsh/ginsh.1.in
@@ -350,6 +350,9 @@ detail here. Please refer to the GiNaC documentation.
 .BI quo( expression ", " expression ", " symbol )
 \- quotient of polynomials
 .br
+.BI rank( matrix )
+\- rank of a matrix
+.br
 .BI rem( expression ", " expression ", " symbol )
 \- remainder of polynomials
 .br

diff --git a/ginsh/ginsh_parser.yy b/ginsh/ginsh_parser.yy
index c045826754f5dc2a4c6f9d0a2f7d0a47aec4852d..ed2fcda6671996e198f0f590f2b2f4e30c65f023 100644 (file)
--- a/ginsh/ginsh_parser.yy
+++ b/ginsh/ginsh_parser.yy
@@ -469,6 +469,12 @@ static ex f_quo(const exprseq &e)
        return quo(e[0], e[1], e[2]);
 }
 
+static ex f_rank(const exprseq &e)
+{
+       CHECK_ARG(0, matrix, rank);
+       return ex_to<matrix>(e[0]).rank();
+}
+
 static ex f_rem(const exprseq &e)
 {
        return rem(e[0], e[1], e[2]);
@@ -581,6 +587,7 @@ static const fcn_init builtin_fcns[] = {
        {"print_csrc", f_dummy, 0},  // for Tab-completion
        {"print_latex", f_dummy, 0}, // for Tab-completion
        {"quo", f_quo, 3},
+       {"rank", f_rank, 1},
        {"rem", f_rem, 3},
        {"series", f_series, 3},
        {"sprem", f_sprem, 3},