author Richard Kreckel Tue, 15 Jun 2004 01:48:54 +0000 (01:48 +0000) committer Richard Kreckel Tue, 15 Jun 2004 01:48:54 +0000 (01:48 +0000)
- Fixes to matrix elimination routines, needed for rank computations.

 check/exam_matrices.cpp patch | blob | history doc/tutorial/ginac.texi patch | blob | history ginac/matrix.cpp patch | blob | history ginac/matrix.h patch | blob | history ginsh/ginsh.1.in patch | blob | history ginsh/ginsh_parser.yy patch | blob | history

index 67c52be..449b6d3 100644 (file)
@@ -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) {
index 3f98183..b58d631 100644 (file)
@@ -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
index 2c6e5bb..e24f644 100644 (file)
@@ -764,7 +764,7 @@ ex matrix::determinant(unsigned algo) const
else
return m.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();
index ee30601..29af2a0 100644 (file)
@@ -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. */
index 01ecbb5..1662ce2 100644 (file)
@@ -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
index c045826..ed2fcda 100644 (file)
@@ -469,6 +469,12 @@ static ex f_quo(const exprseq &e)
return quo(e, e, e);
}

+static ex f_rank(const exprseq &e)
+{
+       CHECK_ARG(0, matrix, rank);
+       return ex_to<matrix>(e).rank();
+}
+
static ex f_rem(const exprseq &e)
{
return rem(e, e, e);
@@ -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},