the destructor, copy constructor, and assignment operator (which were the

[ginac.git] / ginac / matrix.cpp

diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index d6091aa4c1b568c18cdcde9f0944160724e4df1e..0b87ed5d3ba843511e144bc3c316bd99dd0a4d50 100644 (file)
--- 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
@@ -54,28 +54,6 @@ matrix::matrix() : inherited(TINFO_matrix), row(1), col(1)
        m.push_back(_ex0());
 }
 
-matrix::~matrix()
-{
-       debugmsg("matrix destructor",LOGLEVEL_DESTRUCT);
-       destroy(false);
-}
-
-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(true);
-               copy(other);
-       }
-       return *this;
-}
-
 // protected
 
 void matrix::copy(const matrix & other)
@@ -647,7 +625,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 +651,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();
+                                       }
                                }
                        }
                }
@@ -944,11 +925,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)
@@ -1186,7 +1170,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 +1198,6 @@ ex lst_to_matrix(const ex &l)
        return m;
 }
 
-//////////
-// global constants
-//////////
-
-const matrix some_matrix;
-const std::type_info & typeid_matrix = typeid(some_matrix);
-
 #ifndef NO_NAMESPACE_GINAC
 } // namespace GiNaC
 #endif // ndef NO_NAMESPACE_GINAC