- Made determinant_algo (in flags.h) really work.
[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-2000 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 (int size=3; size<20; ++size) {
33         matrix A(size,size);
34         // populate one element in each row:
35         for (int 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 (int 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 (int size=3; size<8; ++size) {
60         matrix A(size,size);
61         for (int r=0; r<size-1; ++r) {
62             // populate one or two elements in each row:
63             for (int 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 (int 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 (int size=3; size<7; ++size) {
94         matrix A(size,size);
95         for (int co=0; co<size-1; ++co) {
96             // populate one or two elements in each row:
97             for (int 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 lines
107         // to guarantee that the determinant is zero:
108         for (int 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 (int size=2; size<6; ++size) {
128         matrix A(size,size);
129         for (int co=0; co<size; ++co) {
130             for (int ro=0; ro<size; ++ro) {
131                 // populate some elements
132                 ex elem = 0;
133                 if (rand()%(size-1) == 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_bareiss = A.determinant(determinant_algo::bareiss);
141         if ((det_gauss-det_laplace).normal() != 0 ||
142             (det_bareiss-det_laplace).normal() != 0) {
143             clog << "Determinant of " << size << "x" << size << " matrix "
144                  << endl << A << endl
145                  << "is inconsistent between different algorithms:" << endl
146                  << "Gauss elimination:   " << det_gauss << endl
147                  << "Minor elimination:   " << det_laplace << endl
148                  << "Fraction-free elim.: " << det_bareiss << endl;
149             ++result;
150         }
151     }
152     
153     return result;
154 }
155
156 unsigned check_matrices(void)
157 {
158     unsigned result = 0;
159     
160     cout << "checking symbolic matrix manipulations" << flush;
161     clog << "---------symbolic matrix manipulations:" << endl;
162     
163     result += integdom_matrix_determinants();  cout << '.' << flush;
164     result += rational_matrix_determinants();  cout << '.' << flush;
165     result += funny_matrix_determinants();  cout << '.' << flush;
166     result += compare_matrix_determinants();  cout << '.' << flush;
167     
168     if (!result) {
169         cout << " passed " << endl;
170         clog << "(no output)" << endl;
171     } else {
172         cout << " failed " << endl;
173     }
174     
175     return result;
176 }