1 /** @file check_matrices.cpp
3 * Here we test manipulations on GiNaC's symbolic matrices. */
6 * GiNaC Copyright (C) 1999-2002 Johannes Gutenberg University Mainz, Germany
8 * This program is free software; you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation; either version 2 of the License, or
11 * (at your option) any later version.
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
18 * You should have received a copy of the GNU General Public License
19 * along with this program; if not, write to the Free Software
20 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25 /* determinants of some sparse symbolic matrices with coefficients in
26 * an integral domain. */
27 static unsigned integdom_matrix_determinants(void)
32 for (unsigned size=3; size<20; ++size) {
34 // populate one element in each row:
35 for (unsigned r=0; r<size-1; ++r)
36 A.set(r,unsigned(rand()%size),dense_univariate_poly(a,5));
37 // set the last row to a linear combination of two other lines
38 // to guarantee that the determinant is zero:
39 for (unsigned c=0; c<size; ++c)
40 A.set(size-1,c,A(0,c)-A(size-2,c));
41 if (!A.determinant().is_zero()) {
42 clog << "Determinant of " << size << "x" << size << " matrix "
44 << "was not found to vanish!" << endl;
52 /* determinants of some symbolic matrices with multivariate rational function
54 static unsigned rational_matrix_determinants(void)
57 symbol a("a"), b("b"), c("c");
59 for (unsigned size=3; size<8; ++size) {
61 for (unsigned r=0; r<size-1; ++r) {
62 // populate one or two elements in each row:
63 for (unsigned ec=0; ec<2; ++ec) {
64 ex numer = sparse_tree(a, b, c, 1+rand()%4, false, false, false);
67 denom = sparse_tree(a, b, c, rand()%2, false, false, false);
68 } while (denom.is_zero());
69 A.set(r,unsigned(rand()%size),numer/denom);
72 // set the last row to a linear combination of two other lines
73 // to guarantee that the determinant is zero:
74 for (unsigned co=0; co<size; ++co)
75 A.set(size-1,co,A(0,co)-A(size-2,co));
76 if (!A.determinant().is_zero()) {
77 clog << "Determinant of " << size << "x" << size << " matrix "
79 << "was not found to vanish!" << endl;
87 /* Some quite funny determinants with functions and stuff like that inside. */
88 static unsigned funny_matrix_determinants(void)
91 symbol a("a"), b("b"), c("c");
93 for (unsigned size=3; size<7; ++size) {
95 for (unsigned co=0; co<size-1; ++co) {
96 // populate one or two elements in each row:
97 for (unsigned ec=0; ec<2; ++ec) {
98 ex numer = sparse_tree(a, b, c, 1+rand()%3, true, true, false);
101 denom = sparse_tree(a, b, c, rand()%2, false, true, false);
102 } while (denom.is_zero());
103 A.set(unsigned(rand()%size),co,numer/denom);
106 // set the last column to a linear combination of two other columns
107 // to guarantee that the determinant is zero:
108 for (unsigned ro=0; ro<size; ++ro)
109 A.set(ro,size-1,A(ro,0)-A(ro,size-2));
110 if (!A.determinant().is_zero()) {
111 clog << "Determinant of " << size << "x" << size << " matrix "
113 << "was not found to vanish!" << endl;
121 /* compare results from different determinant algorithms.*/
122 static unsigned compare_matrix_determinants(void)
127 for (unsigned size=2; size<7; ++size) {
129 for (unsigned co=0; co<size; ++co) {
130 for (unsigned ro=0; ro<size; ++ro) {
131 // populate some elements
133 if (rand()%(size/2) == 0)
134 elem = sparse_tree(a, a, a, rand()%3, false, true, false);
138 ex det_gauss = A.determinant(determinant_algo::gauss);
139 ex det_laplace = A.determinant(determinant_algo::laplace);
140 ex det_divfree = A.determinant(determinant_algo::divfree);
141 ex det_bareiss = A.determinant(determinant_algo::bareiss);
142 if ((det_gauss-det_laplace).normal() != 0 ||
143 (det_bareiss-det_laplace).normal() != 0 ||
144 (det_divfree-det_laplace).normal() != 0) {
145 clog << "Determinant of " << size << "x" << size << " matrix "
147 << "is inconsistent between different algorithms:" << endl
148 << "Gauss elimination: " << det_gauss << endl
149 << "Minor elimination: " << det_laplace << endl
150 << "Division-free elim.: " << det_divfree << endl
151 << "Fraction-free elim.: " << det_bareiss << endl;
159 static unsigned symbolic_matrix_inverse(void)
162 symbol a("a"), b("b"), c("c");
164 for (unsigned size=2; size<5; ++size) {
167 for (unsigned co=0; co<size; ++co) {
168 for (unsigned ro=0; ro<size; ++ro) {
169 // populate some elements
171 if (rand()%(size/2) == 0)
172 elem = sparse_tree(a, b, c, rand()%2, false, true, false);
176 } while (A.determinant() == 0);
177 matrix B = A.inverse();
180 for (unsigned ro=0; ro<size; ++ro)
181 for (unsigned co=0; co<size; ++co)
182 if (C(ro,co).normal() != (ro==co?1:0))
185 clog << "Inverse of " << size << "x" << size << " matrix "
187 << "erroneously returned: "
188 << endl << B << endl;
196 unsigned check_matrices(void)
200 cout << "checking symbolic matrix manipulations" << flush;
201 clog << "---------symbolic matrix manipulations:" << endl;
203 result += integdom_matrix_determinants(); cout << '.' << flush;
204 result += rational_matrix_determinants(); cout << '.' << flush;
205 result += funny_matrix_determinants(); cout << '.' << flush;
206 result += compare_matrix_determinants(); cout << '.' << flush;
207 result += symbolic_matrix_inverse(); cout << '.' << flush;
210 cout << " passed " << endl;
211 clog << "(no output)" << endl;
213 cout << " failed " << endl;