X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=check%2Fcheck_matrices.cpp;h=465c244ca1a00f6a0733473d5bd2e872c109e636;hp=c1fa6257599938bc57c015e363a6822bf876a9c2;hb=7e7beee2c946694130a484c923f6af8391867495;hpb=6db7c98ab607db5329f7c94e3524af67fd80143a

diff --git a/check/check_matrices.cpp b/check/check_matrices.cpp
index c1fa6257..465c244c 100644
--- a/check/check_matrices.cpp
+++ b/check/check_matrices.cpp
@@ -3,7 +3,7 @@
  *  Here we test manipulations on GiNaC's symbolic matrices. */
 
 /*
- *  GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2002 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
@@ -26,115 +26,192 @@
  * an integral domain. */
 static unsigned integdom_matrix_determinants(void)
 {
-    unsigned result = 0;
-    symbol a("a");
-    
-    for (int size=3; size<20; ++size) {
-        matrix A(size,size);
-        for (int r=0; r<size-1; ++r) {
-            // populate one element in each row:
-            A.set(r,unsigned(rand()%size),dense_univariate_poly(a,5));
-        }
-        for (int c=0; c<size; ++c) {
-            // set the last line to a linear combination of two other lines
-            // to guarantee that the determinant is zero:
-            A.set(size-1,c,A(0,c)-A(size-2,c));
-        }
-        if (!A.determinant().is_zero()) {
-            clog << "Determinant of " << size << "x" << size << " matrix "
-                 << endl << A << endl
-                 << "was not found to vanish!" << endl;
-            ++result;
-        }
-    }
-    
-    return result;
+	unsigned result = 0;
+	symbol a("a");
+	
+	for (unsigned size=3; size<20; ++size) {
+		matrix A(size,size);
+		// populate one element in each row:
+		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:
+		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 "
+			     << endl << A << endl
+			     << "was not found to vanish!" << endl;
+			++result;
+		}
+	}
+	
+	return result;
 }
 
-/* determinants of some sparse symbolic matrices with multivariate rational
- * function coefficients. */
+/* determinants of some symbolic matrices with multivariate rational function
+ * coefficients. */
 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 (unsigned size=3; size<8; ++size) {
+		matrix A(size,size);
+		for (unsigned r=0; r<size-1; ++r) {
+			// populate one or two elements in each row:
+			for (unsigned ec=0; ec<2; ++ec) {
+				ex numer = sparse_tree(a, b, c, 1+rand()%4, false, false, false);
+				ex denom;
+				do {
+					denom = sparse_tree(a, b, c, rand()%2, false, false, false);
+				} while (denom.is_zero());
+				A.set(r,unsigned(rand()%size),numer/denom);
+			}
+		}
+		// set the last row to a linear combination of two other lines
+		// to guarantee that the determinant is zero:
+		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 "
+			     << endl << A << endl
+			     << "was not found to vanish!" << endl;
+			++result;
+		}
+	}
+	
+	return result;
+}
 
-    for (int size=3; size<8; ++size) {
-        matrix A(size,size);
-        for (int r=0; r<size-1; ++r) {
-            // populate one element in each row:
-            ex numer = sparse_tree(a, b, c, 4, false, false, false);
-            ex denom;
-            do {
-                denom = sparse_tree(a, b, c, 1, false, false, false);
-            } while (denom.is_zero());
-            A.set(r,unsigned(rand()%size),numer/denom);
-        }
-        for (int c=0; c<size; ++c) {
-            // set the last line to a linear combination of two other lines
-            // to guarantee that the determinant is zero:
-            A.set(size-1,c,A(0,c)-A(size-2,c));
-        }
-        if (!A.determinant().is_zero()) {
-            clog << "Determinant of " << size << "x" << size << " matrix "
-                 << endl << A << endl
-                 << "was not found to vanish!" << endl;
-            ++result;
-        }
-    }
-    
-    return result;
+/* Some quite funny determinants with functions and stuff like that inside. */
+static unsigned funny_matrix_determinants(void)
+{
+	unsigned result = 0;
+	symbol a("a"), b("b"), c("c");
+	
+	for (unsigned size=3; size<7; ++size) {
+		matrix A(size,size);
+		for (unsigned co=0; co<size-1; ++co) {
+			// populate one or two elements in each row:
+			for (unsigned ec=0; ec<2; ++ec) {
+				ex numer = sparse_tree(a, b, c, 1+rand()%3, true, true, false);
+				ex denom;
+				do {
+					denom = sparse_tree(a, b, c, rand()%2, false, true, false);
+				} while (denom.is_zero());
+				A.set(unsigned(rand()%size),co,numer/denom);
+			}
+		}
+		// set the last column to a linear combination of two other columns
+		// to guarantee that the determinant is zero:
+		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 "
+			     << endl << A << endl
+			     << "was not found to vanish!" << endl;
+			++result;
+		}
+	}
+	
+	return result;
+}
+
+/* compare results from different determinant algorithms.*/
+static unsigned compare_matrix_determinants(void)
+{
+	unsigned result = 0;
+	symbol a("a");
+	
+	for (unsigned size=2; size<7; ++size) {
+		matrix A(size,size);
+		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, 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 ||
+			(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
+			     << "Division-free elim.: " << det_divfree << endl
+			     << "Fraction-free elim.: " << det_bareiss << endl;
+			++result;
+		}
+	}
+	
+	return result;
 }
 
-/* Some quite wild determinants with functions and stuff like that. */
-static unsigned wild_matrix_determinants(void)
+static unsigned symbolic_matrix_inverse(void)
 {
-    unsigned result = 0;
-    symbol a("a"), b("b"), c("c");
-    
-    for (int size=3; size<6; ++size) {
-        matrix A(size,size);
-        for (int r=0; r<size-1; ++r) {
-            // populate one element in each row:
-            ex numer = sparse_tree(a, b, c, 3, true, true, false);
-            ex denom;
-            do {
-                denom = sparse_tree(a, b, c, 1, false, true, false);
-            } while (denom.is_zero());
-            A.set(r,unsigned(rand()%size),numer/denom);
-        }
-        for (int c=0; c<size; ++c) {
-            // set the last line to a linear combination of two other lines
-            // to guarantee that the determinant is zero:
-            A.set(size-1,c,A(0,c)-A(size-2,c));
-        }
-        if (!A.determinant().is_zero()) {
-            clog << "Determinant of " << size << "x" << size << " matrix "
-                 << endl << A << endl
-                 << "was not found to vanish!" << endl;
-            ++result;
-        }
-    }
-    
-    return result;
+	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;
-    
-    cout << "checking symbolic matrix manipulations" << flush;
-    clog << "---------symbolic matrix manipulations:" << endl;
-    
-    result += integdom_matrix_determinants();  cout << '.' << flush;
-    result += rational_matrix_determinants();  cout << '.' << flush;
-    result += wild_matrix_determinants();  cout << '.' << flush;
-    
-    if (!result) {
-        cout << " passed " << endl;
-        clog << "(no output)" << endl;
-    } else {
-        cout << " failed " << endl;
-    }
-    
-    return result;
+	unsigned result = 0;
+	
+	cout << "checking symbolic matrix manipulations" << flush;
+	clog << "---------symbolic matrix manipulations:" << endl;
+	
+	result += integdom_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;
+		clog << "(no output)" << endl;
+	} else {
+		cout << " failed " << endl;
+	}
+	
+	return result;
 }