- changed function::diff() to be more tolerant by checking first if the

[ginac.git] / ginac / matrix.cpp

diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index 59863ea5e206f0492ce9268f9006d8ded0582922..bacacad8889a4199c7a07576470903f6586277e7 100644 (file)
--- a/ginac/matrix.cpp
+++ b/ginac/matrix.cpp
@@ -1,7 +1,8 @@
 /** @file matrix.cpp
  *
- *  Implementation of symbolic matrices
- *
+ *  Implementation of symbolic matrices */
+
+/*
  *  GiNaC Copyright (C) 1999 Johannes Gutenberg University Mainz, Germany
  *
  *  This program is free software; you can redistribute it and/or modify
@@ -22,7 +23,10 @@
 #include <algorithm>
 #include <stdexcept>
 
-#include "ginac.h"
+#include "matrix.h"
+#include "debugmsg.h"
+
+namespace GiNaC {
 
 //////////
 // default constructor, destructor, copy constructor, assignment operator
@@ -33,7 +37,7 @@
 
 /** Default ctor.  Initializes to 1 x 1-dimensional zero-matrix. */
 matrix::matrix()
-    : basic(TINFO_MATRIX), row(1), col(1)
+    : basic(TINFO_matrix), row(1), col(1)
 {
     debugmsg("matrix default constructor",LOGLEVEL_CONSTRUCT);
     m.push_back(exZERO());
@@ -86,7 +90,7 @@ void matrix::destroy(bool call_parent)
  *  @param r number of rows
  *  @param c number of cols */
 matrix::matrix(int r, int c)
-    : basic(TINFO_MATRIX), row(r), col(c)
+    : basic(TINFO_matrix), row(r), col(c)
 {
     debugmsg("matrix constructor from int,int",LOGLEVEL_CONSTRUCT);
     m.resize(r*c, exZERO());
@@ -96,7 +100,7 @@ matrix::matrix(int r, int c)
 
 /** Ctor from representation, for internal use only. */
 matrix::matrix(int r, int c, vector<ex> const & m2)
-    : basic(TINFO_MATRIX), row(r), col(c), m(m2)
+    : basic(TINFO_matrix), row(r), col(c), m(m2)
 {
     debugmsg("matrix constructor from int,int,vector<ex>",LOGLEVEL_CONSTRUCT);
 }
@@ -139,7 +143,7 @@ ex matrix::expand(unsigned options) const
  *  itself or one of the elements 'has' it. */
 bool matrix::has(ex const & other) const
 {
-    ASSERT(other.bp!=0);
+    GINAC_ASSERT(other.bp!=0);
     
     // tautology: it is the expression itself
     if (is_equal(*other.bp)) return true;
@@ -209,7 +213,7 @@ ex matrix::evalf(int level) const
 
 int matrix::compare_same_type(basic const & other) const
 {
-    ASSERT(is_exactly_of_type(other, matrix));
+    GINAC_ASSERT(is_exactly_of_type(other, matrix));
     matrix const & o=static_cast<matrix &>(const_cast<basic &>(other));
     
     // compare number of rows
@@ -346,7 +350,7 @@ matrix matrix::transpose(void) const
  * called internally by matrix::determinant(). */
 ex determinant_numeric(const matrix & M)
 {
-    ASSERT(M.rows()==M.cols());  // cannot happen, just in case...
+    GINAC_ASSERT(M.rows()==M.cols());  // cannot happen, just in case...
     matrix tmp(M);
     ex det=exONE();
     ex piv;
@@ -395,7 +399,7 @@ int permutation_sign(vector<T> s)
  *  routine is only called internally by matrix::determinant(). */
 ex determinant_symbolic_perm(const matrix & M)
 {
-    ASSERT(M.rows()==M.cols());  // cannot happen, just in case...
+    GINAC_ASSERT(M.rows()==M.cols());  // cannot happen, just in case...
     
     if (M.rows()==1) {  // speed things up
         return M(0,0);
@@ -420,7 +424,7 @@ ex determinant_symbolic_perm(const matrix & M)
  *  called internally by matrix::determinant(). */
 ex determinant_symbolic_minor(const matrix & M)
 {
-    ASSERT(M.rows()==M.cols());  // cannot happen, just in case...
+    GINAC_ASSERT(M.rows()==M.cols());  // cannot happen, just in case...
     
     if (M.rows()==1) {  // end of recursion
         return M(0,0);
@@ -464,7 +468,7 @@ ex determinant_symbolic_minor(const matrix & M)
  *  that are very hard to canonicalize. */
 /*ex determinant_symbolic_leverrier(const matrix & M)
  *{
- *    ASSERT(M.rows()==M.cols());  // cannot happen, just in case...
+ *    GINAC_ASSERT(M.rows()==M.cols());  // cannot happen, just in case...
  *    
  *    matrix B(M);
  *    matrix I(M.row, M.col);
@@ -686,7 +690,7 @@ matrix matrix::fraction_free_elim(matrix const & vars,
     }
     */
     
-#ifdef DOASSERT
+#ifdef DO_GINAC_ASSERT
     // test if we really have an upper echelon matrix
     int zero_in_last_row=-1;
     for (int r=1; r<=m; ++r) {
@@ -698,10 +702,10 @@ matrix matrix::fraction_free_elim(matrix const & vars,
                 break;
             }
         }
-        ASSERT((zero_in_this_row>zero_in_last_row)||(zero_in_this_row=n));
+        GINAC_ASSERT((zero_in_this_row>zero_in_last_row)||(zero_in_this_row=n));
         zero_in_last_row=zero_in_this_row;
     }
-#endif // def DOASSERT
+#endif // def DO_GINAC_ASSERT
     
     // assemble solution
     matrix sol(n,1);
@@ -743,7 +747,7 @@ matrix matrix::fraction_free_elim(matrix const & vars,
     }
     */
     
-#ifdef DOASSERT
+#ifdef DO_GINAC_ASSERT
     // test solution with echelon matrix
     for (int r=1; r<=m; ++r) {
         ex e=0;
@@ -755,7 +759,7 @@ matrix matrix::fraction_free_elim(matrix const & vars,
             cout << "b.ffe_get(" << r<<",1)=" << b.ffe_get(r,1) << endl;
             cout << "diff=" << (e-b.ffe_get(r,1)).normal() << endl;
         }
-        ASSERT((e-b.ffe_get(r,1)).normal().is_zero());
+        GINAC_ASSERT((e-b.ffe_get(r,1)).normal().is_zero());
     }
 
     // test solution with original matrix
@@ -778,9 +782,9 @@ matrix matrix::fraction_free_elim(matrix const & vars,
             ex xxx=e-rhs.ffe_get(r,1);
             cerr << "xxx=" << xxx << endl << endl;
         }
-        ASSERT((e-rhs.ffe_get(r,1)).normal().is_zero());
+        GINAC_ASSERT((e-rhs.ffe_get(r,1)).normal().is_zero());
     }
-#endif // def DOASSERT
+#endif // def DO_GINAC_ASSERT
     
     return sol;
 }   
@@ -867,3 +871,5 @@ int matrix::pivot(int ro)
 
 const matrix some_matrix;
 type_info const & typeid_matrix=typeid(some_matrix);
+
+} // namespace GiNaC