X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fmatrix.cpp;h=0c32a76796a176be71c55797ec6e438a42abcc92;hp=470991b7bf12d4354afab97ac3226e37514fac92;hb=4e3a4ac2bcb0837611ea31bc8fc05d84a20c33ac;hpb=703c6cebb5d3d395437e73e6935f3691aed68e0a

diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index 470991b7..0c32a767 100644
--- a/ginac/matrix.cpp
+++ b/ginac/matrix.cpp
@@ -3,7 +3,7 @@
  *  Implementation of symbolic matrices */
 
 /*
- *  GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany
  *
  *  This program is free software; you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
@@ -34,50 +34,26 @@
 #include "symbol.h"
 #include "normal.h"
 
-#ifndef NO_NAMESPACE_GINAC
 namespace GiNaC {
-#endif // ndef NO_NAMESPACE_GINAC
 
 GINAC_IMPLEMENT_REGISTERED_CLASS(matrix, basic)
 
 //////////
-// default constructor, destructor, copy constructor, assignment operator
-// and helpers:
+// default ctor, dtor, copy ctor, assignment operator and helpers:
 //////////
 
 // public
 
 /** Default ctor.  Initializes to 1 x 1-dimensional zero-matrix. */
-matrix::matrix()
-	: inherited(TINFO_matrix), row(1), col(1)
+matrix::matrix() : inherited(TINFO_matrix), row(1), col(1)
 {
-	debugmsg("matrix default constructor",LOGLEVEL_CONSTRUCT);
+	debugmsg("matrix default ctor",LOGLEVEL_CONSTRUCT);
 	m.push_back(_ex0());
 }
 
-matrix::~matrix()
-{
-	debugmsg("matrix destructor",LOGLEVEL_DESTRUCT);
-}
-
-matrix::matrix(const matrix & other)
-{
-	debugmsg("matrix copy constructor",LOGLEVEL_CONSTRUCT);
-	copy(other);
-}
-
-const matrix & matrix::operator=(const matrix & other)
-{
-	debugmsg("matrix operator=",LOGLEVEL_ASSIGNMENT);
-	if (this != &other) {
-		destroy(1);
-		copy(other);
-	}
-	return *this;
-}
-
 // protected
 
+/** For use by copy ctor and assignment operator. */
 void matrix::copy(const matrix & other)
 {
 	inherited::copy(other);
@@ -92,7 +68,7 @@ void matrix::destroy(bool call_parent)
 }
 
 //////////
-// other constructors
+// other ctors
 //////////
 
 // public
@@ -102,9 +78,9 @@ void matrix::destroy(bool call_parent)
  *  @param r number of rows
  *  @param c number of cols */
 matrix::matrix(unsigned r, unsigned c)
-	: inherited(TINFO_matrix), row(r), col(c)
+  : inherited(TINFO_matrix), row(r), col(c)
 {
-	debugmsg("matrix constructor from unsigned,unsigned",LOGLEVEL_CONSTRUCT);
+	debugmsg("matrix ctor from unsigned,unsigned",LOGLEVEL_CONSTRUCT);
 	m.resize(r*c, _ex0());
 }
 
@@ -112,9 +88,9 @@ matrix::matrix(unsigned r, unsigned c)
 
 /** Ctor from representation, for internal use only. */
 matrix::matrix(unsigned r, unsigned c, const exvector & m2)
-	: inherited(TINFO_matrix), row(r), col(c), m(m2)
+  : inherited(TINFO_matrix), row(r), col(c), m(m2)
 {
-	debugmsg("matrix constructor from unsigned,unsigned,exvector",LOGLEVEL_CONSTRUCT);
+	debugmsg("matrix ctor from unsigned,unsigned,exvector",LOGLEVEL_CONSTRUCT);
 }
 
 //////////
@@ -124,7 +100,7 @@ matrix::matrix(unsigned r, unsigned c, const exvector & m2)
 /** Construct object from archive_node. */
 matrix::matrix(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
 {
-	debugmsg("matrix constructor from archive_node", LOGLEVEL_CONSTRUCT);
+	debugmsg("matrix ctor from archive_node", LOGLEVEL_CONSTRUCT);
 	if (!(n.find_unsigned("row", row)) || !(n.find_unsigned("col", col)))
 		throw (std::runtime_error("unknown matrix dimensions in archive"));
 	m.reserve(row * col);
@@ -162,12 +138,6 @@ void matrix::archive(archive_node &n) const
 
 // public
 
-basic * matrix::duplicate() const
-{
-	debugmsg("matrix duplicate",LOGLEVEL_DUPLICATE);
-	return new matrix(*this);
-}
-
 void matrix::print(std::ostream & os, unsigned upper_precedence) const
 {
 	debugmsg("matrix print",LOGLEVEL_PRINT);
@@ -393,7 +363,7 @@ matrix matrix::mul(const matrix & other) const
  *  @exception range_error (index out of range) */
 const ex & matrix::operator() (unsigned ro, unsigned co) const
 {
-	if (ro<0 || ro>=row || co<0 || co>=col)
+	if (ro>=row || co>=col)
 		throw (std::range_error("matrix::operator(): index out of range"));
 
 	return m[ro*col+co];
@@ -405,9 +375,9 @@ const ex & matrix::operator() (unsigned ro, unsigned co) const
  *  @exception range_error (index out of range) */
 matrix & matrix::set(unsigned ro, unsigned co, ex value)
 {
-	if (ro<0 || ro>=row || co<0 || co>=col)
+	if (ro>=row || co>=col)
 		throw (std::range_error("matrix::set(): index out of range"));
-	
+    
 	ensure_if_modifiable();
 	m[ro*col+co] = value;
 	return *this;
@@ -647,7 +617,7 @@ matrix matrix::inverse(void) const
 		throw (std::logic_error("matrix::inverse(): matrix not square"));
 	
 	// NOTE: the Gauss-Jordan elimination used here can in principle be
-	// replaced this by two clever calls to gauss_elimination() and some to
+	// replaced by two clever calls to gauss_elimination() and some to
 	// transpose().  Wouldn't be more efficient (maybe less?), just more
 	// orthogonal.
 	matrix tmp(row,col);
@@ -673,14 +643,17 @@ matrix matrix::inverse(void) const
 		}
 		for (unsigned r2=0; r2<row; ++r2) {
 			if (r2 != r1) {
-				ex a2 = cpy.m[r2*col+r1];
-				for (unsigned c=0; c<col; ++c) {
-					cpy.m[r2*col+c] -= a2 * cpy.m[r1*col+c];
-					if (!cpy.m[r2*col+c].info(info_flags::numeric))
-						cpy.m[r2*col+c] = cpy.m[r2*col+c].normal();
-					tmp.m[r2*col+c] -= a2 * tmp.m[r1*col+c];
-					if (!tmp.m[r2*col+c].info(info_flags::numeric))
-						tmp.m[r2*col+c] = tmp.m[r2*col+c].normal();
+				if (!cpy.m[r2*col+r1].is_zero()) {
+					ex a2 = cpy.m[r2*col+r1];
+					// yes, there is something to do in this column
+					for (unsigned c=0; c<col; ++c) {
+						cpy.m[r2*col+c] -= a2 * cpy.m[r1*col+c];
+						if (!cpy.m[r2*col+c].info(info_flags::numeric))
+							cpy.m[r2*col+c] = cpy.m[r2*col+c].normal();
+						tmp.m[r2*col+c] -= a2 * tmp.m[r1*col+c];
+						if (!tmp.m[r2*col+c].info(info_flags::numeric))
+							tmp.m[r2*col+c] = tmp.m[r2*col+c].normal();
+					}
 				}
 			}
 		}
@@ -814,8 +787,8 @@ ex matrix::determinant_minor(void) const
 		return (m[0]*m[3]-m[2]*m[1]).expand();
 	if (n==3)
 		return (m[0]*m[4]*m[8]-m[0]*m[5]*m[7]-
-				m[1]*m[3]*m[8]+m[2]*m[3]*m[7]+
-				m[1]*m[5]*m[6]-m[2]*m[4]*m[6]).expand();
+		        m[1]*m[3]*m[8]+m[2]*m[3]*m[7]+
+		        m[1]*m[5]*m[6]-m[2]*m[4]*m[6]).expand();
 	
 	// This algorithm can best be understood by looking at a naive
 	// implementation of Laplace-expansion, like this one:
@@ -902,7 +875,7 @@ ex matrix::determinant_minor(void) const
 				if (Pkey[fc-1]<fc+c)
 					break;
 			}
-			if (fc<n-c)
+			if (fc<n-c && fc>0)
 				for (unsigned j=fc; j<n-c; ++j)
 					Pkey[j] = Pkey[j-1]+1;
 		} while(fc);
@@ -944,11 +917,14 @@ int matrix::gauss_elimination(const bool det)
 			if (indx > 0)
 				sign = -sign;
 			for (unsigned r2=r0+1; r2<m; ++r2) {
-				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();
+				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)
@@ -1099,15 +1075,15 @@ int matrix::fraction_free_elimination(const bool det)
 			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();
+					              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();
+					              tmp_d.m[r0*n+r1]*tmp_d.m[r2*n+c]).expand();
 					bool check = divide(dividend_n, divisor_n,
-										tmp_n.m[r2*n+c], true);
+					                    tmp_n.m[r2*n+c], true);
 					check &= divide(dividend_d, divisor_d,
-									tmp_d.m[r2*n+c], true);
+					                tmp_d.m[r2*n+c], true);
 					GINAC_ASSERT(check);
 				}
 				// fill up left hand side with zeros
@@ -1186,7 +1162,7 @@ int matrix::pivot(unsigned ro, unsigned co, bool symbolic)
 	// matrix needs pivoting, so swap rows k and ro
 	ensure_if_modifiable();
 	for (unsigned c=0; c<col; ++c)
-		m[k*col+c].swap(m[ro*col+c]);
+		this->m[k*col+c].swap(this->m[ro*col+c]);
 	
 	return k;
 }
@@ -1214,13 +1190,4 @@ ex lst_to_matrix(const ex &l)
 	return m;
 }
 
-//////////
-// global constants
-//////////
-
-const matrix some_matrix;
-const type_info & typeid_matrix=typeid(some_matrix);
-
-#ifndef NO_NAMESPACE_GINAC
 } // namespace GiNaC
-#endif // ndef NO_NAMESPACE_GINAC