Make matrix::solve() work with non-normal zeros.

[ginac.git] / ginac / matrix.cpp

diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index c842885f8103c0a2cc6530b6fc65b6b203197c56..bf7f8347c9a1f7b31645595d408ce602cf77081d 100644 (file)
--- a/ginac/matrix.cpp
+++ b/ginac/matrix.cpp
@@ -3,7 +3,7 @@
  *  Implementation of symbolic matrices */
 
 /*
- *  GiNaC Copyright (C) 1999-2015 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2018 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
@@ -939,10 +939,11 @@ ex matrix::charpoly(const ex & lambda) const
 
 /** Inverse of this matrix.
  *
+ *  @param algo selects the algorithm (one of solve_algo)
  *  @return    the inverted matrix
  *  @exception logic_error (matrix not square)
  *  @exception runtime_error (singular matrix) */
-matrix matrix::inverse() const
+matrix matrix::inverse(unsigned algo) const
 {
        if (row != col)
                throw (std::logic_error("matrix::inverse(): matrix not square"));
@@ -965,7 +966,7 @@ matrix matrix::inverse() const
        
        matrix sol(row,col);
        try {
-               sol = this->solve(vars,identity);
+               sol = this->solve(vars, identity, algo);
        } catch (const std::runtime_error & e) {
            if (e.what()==std::string("matrix::solve(): inconsistent linear system"))
                        throw (std::runtime_error("matrix::inverse(): singular matrix"));
@@ -1053,11 +1054,11 @@ matrix matrix::solve(const matrix & vars,
                unsigned last_assigned_sol = n+1;
                for (int r=m-1; r>=0; --r) {
                        unsigned fnz = 1;    // first non-zero in row
-                       while ((fnz<=n) && (aug.m[r*(n+p)+(fnz-1)].is_zero()))
+                       while ((fnz<=n) && (aug.m[r*(n+p)+(fnz-1)].normal().is_zero()))
                                ++fnz;
                        if (fnz>n) {
                                // row consists only of zeros, corresponding rhs must be 0, too
-                               if (!aug.m[r*(n+p)+n+co].is_zero()) {
+                               if (!aug.m[r*(n+p)+n+co].normal().is_zero()) {
                                        throw (std::runtime_error("matrix::solve(): inconsistent linear system"));
                                }
                        } else {