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

[ginac.git] / ginac / matrix.cpp

diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index 80ca3e459ff42a4002920fd68466054e0398c42a..bf7f8347c9a1f7b31645595d408ce602cf77081d 100644 (file)
--- a/ginac/matrix.cpp
+++ b/ginac/matrix.cpp
@@ -1054,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
                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
                                ++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 {
                                        throw (std::runtime_error("matrix::solve(): inconsistent linear system"));
                                }
                        } else {