465c244ca1a00f6a0733473d5bd2e872c109e636
[ginac.git] / check / check_matrices.cpp
1 /** @file check_matrices.cpp
2  *
3  *  Here we test manipulations on GiNaC's symbolic matrices. */
4
5 /*
6  *  GiNaC Copyright (C) 1999-2002 Johannes Gutenberg University Mainz, Germany
7  *
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.
12  *
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.
17  *
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
21  */
22
23 #include "checks.h"
24
25 /* determinants of some sparse symbolic matrices with coefficients in
26  * an integral domain. */
27 static unsigned integdom_matrix_determinants(void)
28 {
29         unsigned result = 0;
30         symbol a("a");
31         
32         for (unsigned size=3; size<20; ++size) {
33                 matrix A(size,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 "
43                              << endl << A << endl
44                              << "was not found to vanish!" << endl;
45                         ++result;
46                 }
47         }
48         
49         return result;
50 }
51
52 /* determinants of some symbolic matrices with multivariate rational function
53  * coefficients. */
54 static unsigned rational_matrix_determinants(void)
55 {
56         unsigned result = 0;
57         symbol a("a"), b("b"), c("c");
58         
59         for (unsigned size=3; size<8; ++size) {
60                 matrix A(size,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);
65                                 ex denom;
66                                 do {
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);
70                         }
71                 }
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 "
78                              << endl << A << endl
79                              << "was not found to vanish!" << endl;
80                         ++result;
81                 }
82         }
83         
84         return result;
85 }
86
87 /* Some quite funny determinants with functions and stuff like that inside. */
88 static unsigned funny_matrix_determinants(void)
89 {
90         unsigned result = 0;
91         symbol a("a"), b("b"), c("c");
92         
93         for (unsigned size=3; size<7; ++size) {
94                 matrix A(size,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);
99                                 ex denom;
100                                 do {
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);
104                         }
105                 }
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 "
112                              << endl << A << endl
113                              << "was not found to vanish!" << endl;
114                         ++result;
115                 }
116         }
117         
118         return result;
119 }
120
121 /* compare results from different determinant algorithms.*/
122 static unsigned compare_matrix_determinants(void)
123 {
124         unsigned result = 0;
125         symbol a("a");
126         
127         for (unsigned size=2; size<7; ++size) {
128                 matrix A(size,size);
129                 for (unsigned co=0; co<size; ++co) {
130                         for (unsigned ro=0; ro<size; ++ro) {
131                                 // populate some elements
132                                 ex elem = 0;
133                                 if (rand()%(size/2) == 0)
134                                         elem = sparse_tree(a, a, a, rand()%3, false, true, false);
135                                 A.set(ro,co,elem);
136                         }
137                 }
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 "
146                              << endl << A << endl
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;
152                         ++result;
153                 }
154         }
155         
156         return result;
157 }
158
159 static unsigned symbolic_matrix_inverse(void)
160 {
161         unsigned result = 0;
162         symbol a("a"), b("b"), c("c");
163         
164         for (unsigned size=2; size<5; ++size) {
165                 matrix A(size,size);
166                 do {
167                         for (unsigned co=0; co<size; ++co) {
168                                 for (unsigned ro=0; ro<size; ++ro) {
169                                         // populate some elements
170                                         ex elem = 0;
171                                         if (rand()%(size/2) == 0)
172                                                 elem = sparse_tree(a, b, c, rand()%2, false, true, false);
173                                         A.set(ro,co,elem);
174                                 }
175                         }
176                 } while (A.determinant() == 0);
177                 matrix B = A.inverse();
178                 matrix C = A.mul(B);
179                 bool ok = true;
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))
183                                         ok = false;
184                 if (!ok) {
185                         clog << "Inverse of " << size << "x" << size << " matrix "
186                              << endl << A << endl
187                              << "erroneously returned: "
188                              << endl << B << endl;
189                         ++result;
190                 }
191         }
192         
193         return result;
194 }
195
196 unsigned check_matrices(void)
197 {
198         unsigned result = 0;
199         
200         cout << "checking symbolic matrix manipulations" << flush;
201         clog << "---------symbolic matrix manipulations:" << endl;
202         
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;
208         
209         if (!result) {
210                 cout << " passed " << endl;
211                 clog << "(no output)" << endl;
212         } else {
213                 cout << " failed " << endl;
214         }
215         
216         return result;
217 }