3da30af85fcf3870c532e601af2ff3d284d44267
[ginac.git] / check / matrix_checks.cpp
1 /** @file matrix_checks.cpp
2  *
3  *  Here we test manipulations on GiNaC's symbolic matrices. */
4
5 /*
6  *  GiNaC Copyright (C) 1999 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 <stdexcept>
24 #include <ginac/ginac.h>
25
26 #ifndef NO_GINAC_NAMESPACE
27 using namespace GiNaC;
28 #endif // ndef NO_GINAC_NAMESPACE
29
30 static unsigned matrix_determinants(void)
31 {
32     unsigned result = 0;
33     ex det;
34     matrix m1(1,1), m2(2,2), m3(3,3);
35     symbol a("a"), b("b"), c("c");
36     symbol d("d"), e("e"), f("f");
37     symbol g("g"), h("h"), i("i");
38     
39     // check symbolic trivial matrix determinant
40     m1.set(0,0,a);
41     det = m1.determinant();
42     if (det != a) {
43         clog << "determinant of 1x1 matrix " << m1
44              << " erroneously returned " << det << endl;
45         ++result;
46     }
47     
48     // check generic dense symbolic 2x2 matrix determinant
49     m2.set(0,0,a).set(0,1,b);
50     m2.set(1,0,c).set(1,1,d);
51     det = m2.determinant();
52     if (det != (a*d-b*c)) {
53         clog << "determinant of 2x2 matrix " << m2
54              << " erroneously returned " << det << endl;
55         ++result;
56     }
57     
58     // check generic dense symbolic 3x3 matrix determinant
59     m3.set(0,0,a).set(0,1,b).set(0,2,c);
60     m3.set(1,0,d).set(1,1,e).set(1,2,f);
61     m3.set(2,0,g).set(2,1,h).set(2,2,i);
62     det = m3.determinant().expand();
63     if (det != (a*e*i - a*f*h - d*b*i + d*c*h + g*b*f - g*c*e)) {
64         clog << "determinant of 3x3 matrix " << m3
65              << " erroneously returned " << det << endl;
66         ++result;
67     }
68     
69     // check dense numeric 3x3 matrix determinant
70     m3.set(0,0,numeric(0)).set(0,1,numeric(-1)).set(0,2,numeric(3));
71     m3.set(1,0,numeric(3)).set(1,1,numeric(-2)).set(1,2,numeric(2));
72     m3.set(2,0,numeric(3)).set(2,1,numeric(4)).set(2,2,numeric(-2));
73     det = m3.determinant();
74     if (det != 42) {
75         clog << "determinant of 3x3 matrix " << m3
76              << " erroneously returned " << det << endl;
77         ++result;
78     }
79     
80     // check dense symbolic 2x2 matrix determinant
81     m2.set(0,0,a/(a-b)).set(0,1,numeric(1));
82     m2.set(1,0,b/(a-b)).set(1,1,numeric(1));
83     det = m2.determinant(true);
84     if (det != 1) {
85         clog << "determinant of 2x2 matrix " << m2
86              << " erroneously returned " << det << endl;
87         ++result;
88     }
89
90     // check characteristic polynomial
91     m3.set(0,0,a).set(0,1,-2).set(0,2,2);
92     m3.set(1,0,3).set(1,1,a-1).set(1,2,2);
93     m3.set(2,0,3).set(2,1,4).set(2,2,a-3);
94     ex p = m3.charpoly(a);
95     if (p != 0) {
96         clog << "charpoly of 3x3 matrix " << m3
97              << " erroneously returned " << p << endl;
98         ++result;
99     }
100     
101     return result;
102 }
103
104 static unsigned matrix_invert1(void)
105 {
106     matrix m(1,1);
107     symbol a("a");
108
109     m.set(0,0,a);
110     matrix m_i = m.inverse();
111     
112     if (m_i(0,0) != pow(a,-1)) {
113         clog << "inversion of 1x1 matrix " << m
114              << " erroneously returned " << m_i << endl;
115         return 1;
116     }
117     return 0;
118 }
119
120 static unsigned matrix_invert2(void)
121 {
122     matrix m(2,2);
123     symbol a("a"), b("b"), c("c"), d("d");
124     m.set(0,0,a).set(0,1,b);
125     m.set(1,0,c).set(1,1,d);
126     matrix m_i = m.inverse();
127     ex det = m.determinant().expand();
128     
129     if ( (normal(m_i(0,0)*det) != d) ||
130          (normal(m_i(0,1)*det) != -b) ||
131          (normal(m_i(1,0)*det) != -c) ||
132          (normal(m_i(1,1)*det) != a) ) {
133         clog << "inversion of 2x2 matrix " << m
134              << " erroneously returned " << m_i << endl;
135         return 1;
136     }
137     return 0;
138 }
139
140 static unsigned matrix_invert3(void)
141 {
142     matrix m(3,3);
143     symbol a("a"), b("b"), c("c");
144     symbol d("d"), e("e"), f("f");
145     symbol g("g"), h("h"), i("i");
146     m.set(0,0,a).set(0,1,b).set(0,2,c);
147     m.set(1,0,d).set(1,1,e).set(1,2,f);
148     m.set(2,0,g).set(2,1,h).set(2,2,i);
149     matrix m_i = m.inverse();
150     ex det = m.determinant().normal().expand();
151     
152     if ( (normal(m_i(0,0)*det) != (e*i-f*h)) ||
153          (normal(m_i(0,1)*det) != (c*h-b*i)) ||
154          (normal(m_i(0,2)*det) != (b*f-c*e)) ||
155          (normal(m_i(1,0)*det) != (f*g-d*i)) ||
156          (normal(m_i(1,1)*det) != (a*i-c*g)) ||
157          (normal(m_i(1,2)*det) != (c*d-a*f)) ||
158          (normal(m_i(2,0)*det) != (d*h-e*g)) ||
159          (normal(m_i(2,1)*det) != (b*g-a*h)) ||
160          (normal(m_i(2,2)*det) != (a*e-b*d)) ) {
161         clog << "inversion of 3x3 matrix " << m
162              << " erroneously returned " << m_i << endl;
163         return 1;
164     }
165     return 0;
166 }
167
168 static unsigned matrix_misc(void)
169 {
170     unsigned result = 0;
171     matrix m1(2,2);
172     symbol a("a"), b("b"), c("c"), d("d"), e("e"), f("f");
173     m1.set(0,0,a).set(0,1,b);
174     m1.set(1,0,c).set(1,1,d);
175     ex tr = trace(m1);
176     
177     // check a simple trace
178     if (tr.compare(a+d)) {
179         clog << "trace of 2x2 matrix " << m1
180              << " erroneously returned " << tr << endl;
181         ++result;
182     }
183     
184     // and two simple transpositions
185     matrix m2 = transpose(m1);
186     if (m2(0,0) != a || m2(0,1) != c || m2(1,0) != b || m2(1,1) != d) {
187         clog << "transpose of 2x2 matrix " << m1
188              << " erroneously returned " << m2 << endl;
189         ++result;
190     }
191     matrix m3(3,2);
192     m3.set(0,0,a).set(0,1,b);
193     m3.set(1,0,c).set(1,1,d);
194     m3.set(2,0,e).set(2,1,f);
195     if (transpose(transpose(m3)) != m3) {
196         clog << "transposing 3x2 matrix " << m3 << " twice"
197              << " erroneously returned " << transpose(transpose(m3)) << endl;
198         ++result;
199     }
200     
201     // produce a runtime-error by inverting a singular matrix and catch it
202     matrix m4(2,2);
203     matrix m5;
204     bool caught=false;
205     try {
206         m5 = inverse(m4);
207     }
208     catch (std::runtime_error err) {
209         caught=true;
210     }
211     if (!caught) {
212         cerr << "singular 2x2 matrix " << m4
213              << " erroneously inverted to " << m5 << endl;
214         ++result;
215     }
216     
217     return result;
218 }
219
220 unsigned matrix_checks(void)
221 {
222     unsigned result = 0;
223     
224     cout << "checking symbolic matrix manipulations..." << flush;
225     clog << "---------symbolic matrix manipulations:" << endl;
226     
227     result += matrix_determinants();
228     result += matrix_invert1();
229     result += matrix_invert2();
230     result += matrix_invert3();
231     result += matrix_misc();
232     
233     if (! result) {
234         cout << " passed ";
235         clog << "(no output)" << endl;
236     } else {
237         cout << " failed ";
238     }
239     
240     return result;
241 }