location of C++ FAQ Lite has changed
[ginac.git] / check / check_matrices.cpp
1 /** @file check_matrices.cpp
2  *
3  *  Here we test manipulations on GiNaC's symbolic matrices.  They are a
4  *  well-tried resource for cross-checking the underlying symbolic
5  *  manipulations. */
6
7 /*
8  *  GiNaC Copyright (C) 1999-2003 Johannes Gutenberg University Mainz, Germany
9  *
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.
14  *
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.
19  *
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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
23  */
24
25 #include "checks.h"
26
27 /* determinants of some sparse symbolic matrices with coefficients in
28  * an integral domain. */
29 static unsigned integdom_matrix_determinants()
30 {
31         unsigned result = 0;
32         symbol a("a");
33         
34         for (unsigned size=3; size<22; ++size) {
35                 matrix A(size,size);
36                 // populate one element in each row:
37                 for (unsigned r=0; r<size-1; ++r)
38                         A.set(r,unsigned(rand()%size),dense_univariate_poly(a,5));
39                 // set the last row to a linear combination of two other lines
40                 // to guarantee that the determinant is zero:
41                 for (unsigned c=0; c<size; ++c)
42                         A.set(size-1,c,A(0,c)-A(size-2,c));
43                 if (!A.determinant().is_zero()) {
44                         clog << "Determinant of " << size << "x" << size << " matrix "
45                              << endl << A << endl
46                              << "was not found to vanish!" << endl;
47                         ++result;
48                 }
49         }
50         
51         return result;
52 }
53
54 /* determinants of some symbolic matrices with multivariate rational function
55  * coefficients. */
56 static unsigned rational_matrix_determinants()
57 {
58         unsigned result = 0;
59         symbol a("a"), b("b"), c("c");
60         
61         for (unsigned size=3; size<9; ++size) {
62                 matrix A(size,size);
63                 for (unsigned r=0; r<size-1; ++r) {
64                         // populate one or two elements in each row:
65                         for (unsigned ec=0; ec<2; ++ec) {
66                                 ex numer = sparse_tree(a, b, c, 1+rand()%4, false, false, false);
67                                 ex denom;
68                                 do {
69                                         denom = sparse_tree(a, b, c, rand()%2, false, false, false);
70                                 } while (denom.is_zero());
71                                 A.set(r,unsigned(rand()%size),numer/denom);
72                         }
73                 }
74                 // set the last row to a linear combination of two other lines
75                 // to guarantee that the determinant is zero:
76                 for (unsigned co=0; co<size; ++co)
77                         A.set(size-1,co,A(0,co)-A(size-2,co));
78                 if (!A.determinant().is_zero()) {
79                         clog << "Determinant of " << size << "x" << size << " matrix "
80                              << endl << A << endl
81                              << "was not found to vanish!" << endl;
82                         ++result;
83                 }
84         }
85         
86         return result;
87 }
88
89 /* Some quite funny determinants with functions and stuff like that inside. */
90 static unsigned funny_matrix_determinants()
91 {
92         unsigned result = 0;
93         symbol a("a"), b("b"), c("c");
94         
95         for (unsigned size=3; size<8; ++size) {
96                 matrix A(size,size);
97                 for (unsigned co=0; co<size-1; ++co) {
98                         // populate one or two elements in each row:
99                         for (unsigned ec=0; ec<2; ++ec) {
100                                 ex numer = sparse_tree(a, b, c, 1+rand()%3, true, true, false);
101                                 ex denom;
102                                 do {
103                                         denom = sparse_tree(a, b, c, rand()%2, false, true, false);
104                                 } while (denom.is_zero());
105                                 A.set(unsigned(rand()%size),co,numer/denom);
106                         }
107                 }
108                 // set the last column to a linear combination of two other columns
109                 // to guarantee that the determinant is zero:
110                 for (unsigned ro=0; ro<size; ++ro)
111                         A.set(ro,size-1,A(ro,0)-A(ro,size-2));
112                 if (!A.determinant().is_zero()) {
113                         clog << "Determinant of " << size << "x" << size << " matrix "
114                              << endl << A << endl
115                              << "was not found to vanish!" << endl;
116                         ++result;
117                 }
118         }
119         
120         return result;
121 }
122
123 /* compare results from different determinant algorithms.*/
124 static unsigned compare_matrix_determinants()
125 {
126         unsigned result = 0;
127         symbol a("a");
128         
129         for (unsigned size=2; size<8; ++size) {
130                 matrix A(size,size);
131                 for (unsigned co=0; co<size; ++co) {
132                         for (unsigned ro=0; ro<size; ++ro) {
133                                 // populate some elements
134                                 ex elem = 0;
135                                 if (rand()%(size/2) == 0)
136                                         elem = sparse_tree(a, a, a, rand()%3, false, true, false);
137                                 A.set(ro,co,elem);
138                         }
139                 }
140                 ex det_gauss = A.determinant(determinant_algo::gauss);
141                 ex det_laplace = A.determinant(determinant_algo::laplace);
142                 ex det_divfree = A.determinant(determinant_algo::divfree);
143                 ex det_bareiss = A.determinant(determinant_algo::bareiss);
144                 if ((det_gauss-det_laplace).normal() != 0 ||
145                         (det_bareiss-det_laplace).normal() != 0 ||
146                         (det_divfree-det_laplace).normal() != 0) {
147                         clog << "Determinant of " << size << "x" << size << " matrix "
148                              << endl << A << endl
149                              << "is inconsistent between different algorithms:" << endl
150                              << "Gauss elimination:   " << det_gauss << endl
151                              << "Minor elimination:   " << det_laplace << endl
152                              << "Division-free elim.: " << det_divfree << endl
153                              << "Fraction-free elim.: " << det_bareiss << endl;
154                         ++result;
155                 }
156         }
157         
158         return result;
159 }
160
161 static unsigned symbolic_matrix_inverse()
162 {
163         unsigned result = 0;
164         symbol a("a"), b("b"), c("c");
165         
166         for (unsigned size=2; size<6; ++size) {
167                 matrix A(size,size);
168                 do {
169                         for (unsigned co=0; co<size; ++co) {
170                                 for (unsigned ro=0; ro<size; ++ro) {
171                                         // populate some elements
172                                         ex elem = 0;
173                                         if (rand()%(size/2) == 0)
174                                                 elem = sparse_tree(a, b, c, rand()%2, false, true, false);
175                                         A.set(ro,co,elem);
176                                 }
177                         }
178                 } while (A.determinant() == 0);
179                 matrix B = A.inverse();
180                 matrix C = A.mul(B);
181                 bool ok = true;
182                 for (unsigned ro=0; ro<size; ++ro)
183                         for (unsigned co=0; co<size; ++co)
184                                 if (C(ro,co).normal() != (ro==co?1:0))
185                                         ok = false;
186                 if (!ok) {
187                         clog << "Inverse of " << size << "x" << size << " matrix "
188                              << endl << A << endl
189                              << "erroneously returned: "
190                              << endl << B << endl;
191                         ++result;
192                 }
193         }
194         
195         return result;
196 }
197
198 unsigned check_matrices()
199 {
200         unsigned result = 0;
201         
202         cout << "checking symbolic matrix manipulations" << flush;
203         clog << "---------symbolic matrix manipulations:" << endl;
204         
205         result += integdom_matrix_determinants();  cout << '.' << flush;
206         result += rational_matrix_determinants();  cout << '.' << flush;
207         result += funny_matrix_determinants();  cout << '.' << flush;
208         result += compare_matrix_determinants();  cout << '.' << flush;
209         result += symbolic_matrix_inverse();  cout << '.' << flush;
210         
211         if (!result) {
212                 cout << " passed " << endl;
213                 clog << "(no output)" << endl;
214         } else {
215                 cout << " failed " << endl;
216         }
217         
218         return result;
219 }