Finalize 1.7.6 release.

[ginac.git] / check / exam_matrices.cpp

diff --git a/check/exam_matrices.cpp b/check/exam_matrices.cpp
index 97f0c01073841304661ead3c6ac5b292350e83e7..b163386fe4f22b5d3e8af952f16401649637b65d 100644 (file)
--- a/check/exam_matrices.cpp
+++ b/check/exam_matrices.cpp
@@ -1,9 +1,10 @@
+
 /** @file exam_matrices.cpp
  *
  *  Here we examine manipulations on GiNaC's symbolic matrices. */
 
 /*
- *  GiNaC Copyright (C) 1999-2016 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2019 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
@@ -214,6 +215,31 @@ static unsigned matrix_solve2()
        return result;
 }
 
+static unsigned matrix_solve3()
+{
+       unsigned result = 0;
+       symbol x("x");
+       symbol t1("t1"), t2("t2"), t3("t3");
+       matrix A = {
+               {3+6*x, 6*(x+x*x)/(2+3*x), 0},
+               {-(2+7*x+6*x*x)/x, -2-2*x, 0},
+               {-2*(2+3*x)/(1+2*x), -6*x/(1+2*x), 1+4*x}
+       };
+       matrix B = {{0}, {0}, {0}};
+       matrix X = {{t1}, {t2}, {t3}};
+       for (auto algo : vector<int>({
+               solve_algo::gauss, solve_algo::divfree, solve_algo::bareiss, solve_algo::markowitz
+       })) {
+               matrix sol(A.solve(X, B, algo));
+               if (!normal((A*sol - B).evalm()).is_zero_matrix()) {
+                       clog << "Solving " << A << " * " << X << " == " << B << " with algo=" << algo << endl
+                            << "erroneously returned " << sol << endl;
+                       result += 1;
+               }
+       }
+       return result;
+}
+
 static unsigned matrix_evalm()
 {
        unsigned result = 0;
@@ -275,6 +301,35 @@ static unsigned matrix_rank()
        return result;  
 }
 
+unsigned matrix_solve_nonnormal()
+{
+       symbol a("a"), b("b"), c("c"), x("x");
+       // This matrix has a non-normal zero element!
+       matrix mx {{1,0,0},
+                  {0,1/(x+1)-(x-1)/(x*x-1),1},
+                  {0,0,0}};
+       matrix zero {{0}, {0}, {0}};
+       matrix vars {{a}, {b}, {c}};
+       try {
+               matrix sol_gauss   = mx.solve(vars, zero, solve_algo::gauss);
+               matrix sol_divfree = mx.solve(vars, zero, solve_algo::divfree);
+               matrix sol_bareiss = mx.solve(vars, zero, solve_algo::bareiss);
+               if (sol_gauss != sol_divfree  ||  sol_gauss != sol_bareiss) {
+                       clog << "different solutions while solving "
+                            << mx << " * " << vars << " == " << zero << endl
+                            << "gauss:   " << sol_gauss << endl
+                            << "divfree: " << sol_divfree << endl
+                            << "bareiss: " << sol_bareiss << endl;
+                       return 1;
+               }
+       } catch (const exception & e) {
+               clog << "exception thrown while solving "
+                    << mx << " * " << vars << " == " << zero << endl;
+               return 1;
+       }
+       return 0;
+}
+
 static unsigned matrix_misc()
 {
        unsigned result = 0;
@@ -335,8 +390,10 @@ unsigned exam_matrices()
        result += matrix_invert2();  cout << '.' << flush;
        result += matrix_invert3();  cout << '.' << flush;
        result += matrix_solve2();  cout << '.' << flush;
+       result += matrix_solve3();  cout << '.' << flush;
        result += matrix_evalm();  cout << "." << flush;
        result += matrix_rank();  cout << "." << flush;
+       result += matrix_solve_nonnormal();  cout << "." << flush;
        result += matrix_misc();  cout << '.' << flush;
        
        return result;