From f6a16c5f0c04d44e61870b834f6b2e88c98cc8a2 Mon Sep 17 00:00:00 2001 From: Richard Kreckel Date: Sat, 6 May 2000 19:43:12 +0000 Subject: [PATCH] - matrix::determinant_numeric(): deleted. It was duplicate code because we have matrix::gauss_elimination(). --- ginac/matrix.cpp | 69 ++++++++++++++---------------------------------- ginac/matrix.h | 1 - 2 files changed, 20 insertions(+), 50 deletions(-) diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp index c5d056c2..bb1850f7 100644 --- a/ginac/matrix.cpp +++ b/ginac/matrix.cpp @@ -456,6 +456,7 @@ ex matrix::determinant(void) const if (this->row==1) // continuation would be pointless return m[0]; + // Gather some information about the matrix: bool numeric_flag = true; bool normal_flag = false; unsigned sparse_count = 0; // count non-zero elements @@ -469,9 +470,17 @@ ex matrix::determinant(void) const normal_flag = true; } - if (numeric_flag) // purely numeric matrix - return determinant_numeric(); - // Does anybody really know when a matrix is sparse? + // Purely numeric matrix handled by Gauss elimination + if (numeric_flag) { + ex det = 1; + matrix tmp(*this); + int sign = tmp.gauss_elimination(); + for (int d=0; d 0) sign = -sign; for (unsigned r2=r1+1; r2m[r2*col+r1] / this->m[r1*col+r1]; for (unsigned c=r1+1; cm[r2*col+c] -= this->m[r2*col+r1]*this->m[r1*col+c]/this->m[r1*col+r1]; + this->m[r2*col+c] -= piv * this->m[r1*col+c]; for (unsigned c=0; c<=r1; ++c) this->m[r2*col+c] = _ex0(); } @@ -1130,7 +1112,6 @@ int matrix::fraction_free_elimination(bool det) } for (unsigned r1=0; r1" << string(60,'=') << endl; int indx = tmp_n.pivot(r1); if (det && indx==-1) return 0; // FIXME: what to do if det is false? @@ -1140,8 +1121,6 @@ int matrix::fraction_free_elimination(bool det) for (unsigned c=0; c0) { divisor_n = tmp_n.m[(r1-1)*col+(r1-1)].expand(); divisor_d = tmp_d.m[(r1-1)*col+(r1-1)].expand(); @@ -1161,12 +1140,6 @@ int matrix::fraction_free_elimination(bool det) tmp_d.m[r1*col+r1]*tmp_d.m[r2*col+c]).expand(); dividend_d = (tmp_d.m[r2*col+r1]*tmp_d.m[r1*col+c]* tmp_d.m[r1*col+r1]*tmp_d.m[r2*col+c]).expand(); -// cout << "Element " << r2 << ',' << c << endl; -// cout << "dividend_n==" << dividend_n << endl; -// cout << "dividend_d==" << dividend_d << endl; -// cout << " divisor_n==" << divisor_n << endl; -// cout << " divisor_d==" << divisor_d << endl; -// cout << string(20,'-') << endl; bool check = divide(dividend_n, divisor_n, tmp_n.m[r2*col+c],true); check &= divide(dividend_d, divisor_d, @@ -1177,8 +1150,6 @@ int matrix::fraction_free_elimination(bool det) for (unsigned c=0; c<=r1; ++c) tmp_n.m[r2*col+c] = _ex0(); } -// cout << tmp_n << endl; -// cout << tmp_d << endl; } // repopulate *this matrix: it = m.begin(); @@ -1191,7 +1162,7 @@ int matrix::fraction_free_elimination(bool det) } -/** Partial pivoting method. +/** Partial pivoting method for matrix elimination schemes. * Usual pivoting (symbolic==false) returns the index to the element with the * largest absolute value in column ro and swaps the current row with the one * where the element was found. With (symbolic==true) it does the same thing diff --git a/ginac/matrix.h b/ginac/matrix.h index 9742f398..208b6303 100644 --- a/ginac/matrix.h +++ b/ginac/matrix.h @@ -95,7 +95,6 @@ public: matrix solve(const matrix & vars, const matrix & rhs) const; matrix old_solve(const matrix & v) const; // FIXME: may be removed protected: - ex determinant_numeric(void) const; ex determinant_minor(void) const; int gauss_elimination(void); int division_free_elimination(void); -- 2.25.4