- Complete revamp of methods in class matrix. Some redundant (and poor)

[ginac.git] / check / check_matrices.cpp

diff --git a/check/check_matrices.cpp b/check/check_matrices.cpp
index 3e702b8e963c174bdf37a896b91f89d7f5003d18..62b74ba2651ff5d4134dbed97c380b328cb409c3 100644 (file)
--- a/check/check_matrices.cpp
+++ b/check/check_matrices.cpp
@@ -29,14 +29,14 @@ static unsigned integdom_matrix_determinants(void)
     unsigned result = 0;
     symbol a("a");
     
     unsigned result = 0;
     symbol a("a");
     

-    for (int size=3; size<20; ++size) {
+    for (unsigned size=3; size<20; ++size) {

         matrix A(size,size);
         // populate one element in each row:
         matrix A(size,size);
         // populate one element in each row:

-        for (int r=0; r<size-1; ++r)
+        for (unsigned r=0; r<size-1; ++r)

             A.set(r,unsigned(rand()%size),dense_univariate_poly(a,5));
         // set the last row to a linear combination of two other lines
         // to guarantee that the determinant is zero:
             A.set(r,unsigned(rand()%size),dense_univariate_poly(a,5));
         // set the last row to a linear combination of two other lines
         // to guarantee that the determinant is zero:

-        for (int c=0; c<size; ++c)
+        for (unsigned c=0; c<size; ++c)

             A.set(size-1,c,A(0,c)-A(size-2,c));
         if (!A.determinant().is_zero()) {
             clog << "Determinant of " << size << "x" << size << " matrix "
             A.set(size-1,c,A(0,c)-A(size-2,c));
         if (!A.determinant().is_zero()) {
             clog << "Determinant of " << size << "x" << size << " matrix "


@@ -56,11 +56,11 @@ static unsigned rational_matrix_determinants(void)
     unsigned result = 0;
     symbol a("a"), b("b"), c("c");
     
     unsigned result = 0;
     symbol a("a"), b("b"), c("c");
     

-    for (int size=3; size<8; ++size) {
+    for (unsigned size=3; size<8; ++size) {

         matrix A(size,size);
         matrix A(size,size);

-        for (int r=0; r<size-1; ++r) {
+        for (unsigned r=0; r<size-1; ++r) {

             // populate one or two elements in each row:
             // populate one or two elements in each row:

-            for (int ec=0; ec<2; ++ec) {
+            for (unsigned ec=0; ec<2; ++ec) {

                 ex numer = sparse_tree(a, b, c, 1+rand()%4, false, false, false);
                 ex denom;
                 do {
                 ex numer = sparse_tree(a, b, c, 1+rand()%4, false, false, false);
                 ex denom;
                 do {


@@ -71,7 +71,7 @@ static unsigned rational_matrix_determinants(void)
         }
         // set the last row to a linear combination of two other lines
         // to guarantee that the determinant is zero:
         }
         // set the last row to a linear combination of two other lines
         // to guarantee that the determinant is zero:

-        for (int co=0; co<size; ++co)
+        for (unsigned co=0; co<size; ++co)

             A.set(size-1,co,A(0,co)-A(size-2,co));
         if (!A.determinant().is_zero()) {
             clog << "Determinant of " << size << "x" << size << " matrix "
             A.set(size-1,co,A(0,co)-A(size-2,co));
         if (!A.determinant().is_zero()) {
             clog << "Determinant of " << size << "x" << size << " matrix "


@@ -90,11 +90,11 @@ static unsigned funny_matrix_determinants(void)
     unsigned result = 0;
     symbol a("a"), b("b"), c("c");
     
     unsigned result = 0;
     symbol a("a"), b("b"), c("c");
     

-    for (int size=3; size<7; ++size) {
+    for (unsigned size=3; size<7; ++size) {

         matrix A(size,size);
         matrix A(size,size);

-        for (int co=0; co<size-1; ++co) {
+        for (unsigned co=0; co<size-1; ++co) {

             // populate one or two elements in each row:
             // populate one or two elements in each row:

-            for (int ec=0; ec<2; ++ec) {
+            for (unsigned ec=0; ec<2; ++ec) {

                 ex numer = sparse_tree(a, b, c, 1+rand()%3, true, true, false);
                 ex denom;
                 do {
                 ex numer = sparse_tree(a, b, c, 1+rand()%3, true, true, false);
                 ex denom;
                 do {


@@ -103,9 +103,9 @@ static unsigned funny_matrix_determinants(void)
                 A.set(unsigned(rand()%size),co,numer/denom);
             }
         }
                 A.set(unsigned(rand()%size),co,numer/denom);
             }
         }

-        // set the last column to a linear combination of two other lines
+        // set the last column to a linear combination of two other columns

         // to guarantee that the determinant is zero:
         // to guarantee that the determinant is zero:

-        for (int ro=0; ro<size; ++ro)
+        for (unsigned ro=0; ro<size; ++ro)

             A.set(ro,size-1,A(ro,0)-A(ro,size-2));
         if (!A.determinant().is_zero()) {
             clog << "Determinant of " << size << "x" << size << " matrix "
             A.set(ro,size-1,A(ro,0)-A(ro,size-2));
         if (!A.determinant().is_zero()) {
             clog << "Determinant of " << size << "x" << size << " matrix "


@@ -124,27 +124,30 @@ static unsigned compare_matrix_determinants(void)
     unsigned result = 0;
     symbol a("a");
     
     unsigned result = 0;
     symbol a("a");
     

-    for (int size=2; size<6; ++size) {
+    for (unsigned size=2; size<7; ++size) {

         matrix A(size,size);
         matrix A(size,size);

-        for (int co=0; co<size; ++co) {
-            for (int ro=0; ro<size; ++ro) {
+        for (unsigned co=0; co<size; ++co) {
+            for (unsigned ro=0; ro<size; ++ro) {

                 // populate some elements
                 ex elem = 0;
                 // populate some elements
                 ex elem = 0;

-                if (rand()%(size-1) == 0)
+                if (rand()%(size/2) == 0)

                     elem = sparse_tree(a, a, a, rand()%3, false, true, false);
                 A.set(ro,co,elem);
             }
         }
         ex det_gauss = A.determinant(determinant_algo::gauss);
         ex det_laplace = A.determinant(determinant_algo::laplace);
                     elem = sparse_tree(a, a, a, rand()%3, false, true, false);
                 A.set(ro,co,elem);
             }
         }
         ex det_gauss = A.determinant(determinant_algo::gauss);
         ex det_laplace = A.determinant(determinant_algo::laplace);

+        ex det_divfree = A.determinant(determinant_algo::divfree);

         ex det_bareiss = A.determinant(determinant_algo::bareiss);
         if ((det_gauss-det_laplace).normal() != 0 ||
         ex det_bareiss = A.determinant(determinant_algo::bareiss);
         if ((det_gauss-det_laplace).normal() != 0 ||

-            (det_bareiss-det_laplace).normal() != 0) {
+            (det_bareiss-det_laplace).normal() != 0 ||
+            (det_divfree-det_laplace).normal() != 0) {

             clog << "Determinant of " << size << "x" << size << " matrix "
                  << endl << A << endl
                  << "is inconsistent between different algorithms:" << endl
                  << "Gauss elimination:   " << det_gauss << endl
                  << "Minor elimination:   " << det_laplace << endl
             clog << "Determinant of " << size << "x" << size << " matrix "
                  << endl << A << endl
                  << "is inconsistent between different algorithms:" << endl
                  << "Gauss elimination:   " << det_gauss << endl
                  << "Minor elimination:   " << det_laplace << endl

+                 << "Division-free elim.: " << det_divfree << endl

                  << "Fraction-free elim.: " << det_bareiss << endl;
             ++result;
         }
                  << "Fraction-free elim.: " << det_bareiss << endl;
             ++result;
         }


@@ -153,6 +156,43 @@ static unsigned compare_matrix_determinants(void)
     return result;
 }
 
     return result;
 }
 

+static unsigned symbolic_matrix_inverse(void)
+{
+    unsigned result = 0;
+    symbol a("a"), b("b"), c("c");
+    
+    for (unsigned size=2; size<5; ++size) {
+        matrix A(size,size);
+        do {
+            for (unsigned co=0; co<size; ++co) {
+                for (unsigned ro=0; ro<size; ++ro) {
+                    // populate some elements
+                    ex elem = 0;
+                    if (rand()%(size/2) == 0)
+                        elem = sparse_tree(a, b, c, rand()%2, false, true, false);
+                    A.set(ro,co,elem);
+                }
+            }
+        } while (A.determinant() == 0);
+        matrix B = A.inverse();
+        matrix C = A.mul(B);
+        bool ok = true;
+        for (unsigned ro=0; ro<size; ++ro)
+            for (unsigned co=0; co<size; ++co)
+                if (C(ro,co).normal() != (ro==co?1:0))
+                    ok = false;
+        if (!ok) {
+            clog << "Inverse of " << size << "x" << size << " matrix "
+                 << endl << A << endl
+                 << "erroneously returned: "
+                 << endl << B << endl;
+            ++result;
+        }
+    }
+    
+    return result;
+}
+

 unsigned check_matrices(void)
 {
     unsigned result = 0;
 unsigned check_matrices(void)
 {
     unsigned result = 0;


@@ -164,6 +204,7 @@ unsigned check_matrices(void)
     result += rational_matrix_determinants();  cout << '.' << flush;
     result += funny_matrix_determinants();  cout << '.' << flush;
     result += compare_matrix_determinants();  cout << '.' << flush;
     result += rational_matrix_determinants();  cout << '.' << flush;
     result += funny_matrix_determinants();  cout << '.' << flush;
     result += compare_matrix_determinants();  cout << '.' << flush;

+    result += symbolic_matrix_inverse();  cout << '.' << flush;

     
     if (!result) {
         cout << " passed " << endl;
     
     if (!result) {
         cout << " passed " << endl;