1 /** @file check_matrices.cpp
3 * Here we test manipulations on GiNaC's symbolic matrices. They are a
4 * well-tried resource for cross-checking the underlying symbolic
8 * GiNaC Copyright (C) 1999-2011 Johannes Gutenberg University Mainz, Germany
10 * This program is free software; you can redistribute it and/or modify
11 * it under the terms of the GNU General Public License as published by
12 * the Free Software Foundation; either version 2 of the License, or
13 * (at your option) any later version.
15 * This program is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 * GNU General Public License for more details.
20 * You should have received a copy of the GNU General Public License
21 * along with this program; if not, write to the Free Software
22 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
26 using namespace GiNaC;
28 #include <cstdlib> // for rand(), RAND_MAX
33 sparse_tree(const symbol & x, const symbol & y, const symbol & z,
34 int level, bool trig = false, bool rational = true,
35 bool complex = false);
37 dense_univariate_poly(const symbol & x, unsigned degree);
39 /* determinants of some sparse symbolic matrices with coefficients in
40 * an integral domain. */
41 static unsigned integdom_matrix_determinants()
46 for (unsigned size=3; size<22; ++size) {
48 // populate one element in each row:
49 for (unsigned r=0; r<size-1; ++r)
50 A.set(r,unsigned(rand()%size),dense_univariate_poly(a,5));
51 // set the last row to a linear combination of two other lines
52 // to guarantee that the determinant is zero:
53 for (unsigned c=0; c<size; ++c)
54 A.set(size-1,c,A(0,c)-A(size-2,c));
55 if (!A.determinant().is_zero()) {
56 clog << "Determinant of " << size << "x" << size << " matrix "
58 << "was not found to vanish!" << endl;
66 /* determinants of some symbolic matrices with multivariate rational function
68 static unsigned rational_matrix_determinants()
71 symbol a("a"), b("b"), c("c");
73 for (unsigned size=3; size<9; ++size) {
75 for (unsigned r=0; r<size-1; ++r) {
76 // populate one or two elements in each row:
77 for (unsigned ec=0; ec<2; ++ec) {
78 ex numer = sparse_tree(a, b, c, 1+rand()%4, false, false, false);
81 denom = sparse_tree(a, b, c, rand()%2, false, false, false);
82 } while (denom.is_zero());
83 A.set(r,unsigned(rand()%size),numer/denom);
86 // set the last row to a linear combination of two other lines
87 // to guarantee that the determinant is zero:
88 for (unsigned co=0; co<size; ++co)
89 A.set(size-1,co,A(0,co)-A(size-2,co));
90 if (!A.determinant().is_zero()) {
91 clog << "Determinant of " << size << "x" << size << " matrix "
93 << "was not found to vanish!" << endl;
101 /* Some quite funny determinants with functions and stuff like that inside. */
102 static unsigned funny_matrix_determinants()
105 symbol a("a"), b("b"), c("c");
107 for (unsigned size=3; size<8; ++size) {
109 for (unsigned co=0; co<size-1; ++co) {
110 // populate one or two elements in each row:
111 for (unsigned ec=0; ec<2; ++ec) {
112 ex numer = sparse_tree(a, b, c, 1+rand()%3, true, true, false);
115 denom = sparse_tree(a, b, c, rand()%2, false, true, false);
116 } while (denom.is_zero());
117 A.set(unsigned(rand()%size),co,numer/denom);
120 // set the last column to a linear combination of two other columns
121 // to guarantee that the determinant is zero:
122 for (unsigned ro=0; ro<size; ++ro)
123 A.set(ro,size-1,A(ro,0)-A(ro,size-2));
124 if (!A.determinant().is_zero()) {
125 clog << "Determinant of " << size << "x" << size << " matrix "
127 << "was not found to vanish!" << endl;
135 /* compare results from different determinant algorithms.*/
136 static unsigned compare_matrix_determinants()
141 for (unsigned size=2; size<8; ++size) {
143 for (unsigned co=0; co<size; ++co) {
144 for (unsigned ro=0; ro<size; ++ro) {
145 // populate some elements
147 if (rand()%(size/2) == 0)
148 elem = sparse_tree(a, a, a, rand()%3, false, true, false);
152 ex det_gauss = A.determinant(determinant_algo::gauss);
153 ex det_laplace = A.determinant(determinant_algo::laplace);
154 ex det_divfree = A.determinant(determinant_algo::divfree);
155 ex det_bareiss = A.determinant(determinant_algo::bareiss);
156 if ((det_gauss-det_laplace).normal() != 0 ||
157 (det_bareiss-det_laplace).normal() != 0 ||
158 (det_divfree-det_laplace).normal() != 0) {
159 clog << "Determinant of " << size << "x" << size << " matrix "
161 << "is inconsistent between different algorithms:" << endl
162 << "Gauss elimination: " << det_gauss << endl
163 << "Minor elimination: " << det_laplace << endl
164 << "Division-free elim.: " << det_divfree << endl
165 << "Fraction-free elim.: " << det_bareiss << endl;
173 static unsigned symbolic_matrix_inverse()
176 symbol a("a"), b("b"), c("c");
178 for (unsigned size=2; size<6; ++size) {
181 for (unsigned co=0; co<size; ++co) {
182 for (unsigned ro=0; ro<size; ++ro) {
183 // populate some elements
185 if (rand()%(size/2) == 0)
186 elem = sparse_tree(a, b, c, rand()%2, false, true, false);
190 } while (A.determinant() == 0);
191 matrix B = A.inverse();
194 for (unsigned ro=0; ro<size; ++ro)
195 for (unsigned co=0; co<size; ++co)
196 if (C(ro,co).normal() != (ro==co?1:0))
199 clog << "Inverse of " << size << "x" << size << " matrix "
201 << "erroneously returned: "
202 << endl << B << endl;
210 unsigned check_matrices()
214 cout << "checking symbolic matrix manipulations" << flush;
216 result += integdom_matrix_determinants(); cout << '.' << flush;
217 result += rational_matrix_determinants(); cout << '.' << flush;
218 result += funny_matrix_determinants(); cout << '.' << flush;
219 result += compare_matrix_determinants(); cout << '.' << flush;
220 result += symbolic_matrix_inverse(); cout << '.' << flush;
225 int main(int argc, char** argv)
227 return check_matrices();