- enforced GiNaC coding standards :-)
[ginac.git] / check / matrix_checks.cpp
1 /** @file matrix_checks.cpp
2  *
3  *  Here we test manipulations on GiNaC's symbolic matrices.
4  *
5  *  GiNaC Copyright (C) 1999 Johannes Gutenberg University Mainz, Germany
6  *
7  *  This program is free software; you can redistribute it and/or modify
8  *  it under the terms of the GNU General Public License as published by
9  *  the Free Software Foundation; either version 2 of the License, or
10  *  (at your option) any later version.
11  *
12  *  This program is distributed in the hope that it will be useful,
13  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
14  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15  *  GNU General Public License for more details.
16  *
17  *  You should have received a copy of the GNU General Public License
18  *  along with this program; if not, write to the Free Software
19  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
20  */
21
22 #include <stdexcept>
23 #include <ginac/ginac.h>
24
25 static unsigned matrix_determinants(void)
26 {
27     unsigned result = 0;
28     ex det;
29     matrix m1(1,1), m2(2,2), m3(3,3);
30     symbol a("a"), b("b"), c("c");
31     symbol d("d"), e("e"), f("f");
32     symbol g("g"), h("h"), i("i");
33     
34     // check symbolic trivial matrix determinant
35     m1.set(0,0,a);
36     det = m1.determinant();
37     if (det != a) {
38         clog << "determinant of 1x1 matrix " << m1
39              << " erroneously returned " << det << endl;
40         ++result;
41     }
42     
43     // check generic dense symbolic 2x2 matrix determinant
44     m2.set(0,0,a).set(0,1,b);
45     m2.set(1,0,c).set(1,1,d);
46     det = m2.determinant();
47     if (det != (a*d-b*c)) {
48         clog << "determinant of 2x2 matrix " << m2
49              << " erroneously returned " << det << endl;
50         ++result;
51     }
52     
53     // check generic dense symbolic 3x3 matrix determinant
54     m3.set(0,0,a).set(0,1,b).set(0,2,c);
55     m3.set(1,0,d).set(1,1,e).set(1,2,f);
56     m3.set(2,0,g).set(2,1,h).set(2,2,i);
57     det = m3.determinant().expand();
58     if (det != (a*e*i - a*f*h - d*b*i + d*c*h + g*b*f - g*c*e)) {
59         clog << "determinant of 3x3 matrix " << m3
60              << " erroneously returned " << det << endl;
61         ++result;
62     }
63     
64     // check dense numeric 3x3 matrix determinant
65     m3.set(0,0,numeric(0)).set(0,1,numeric(-1)).set(0,2,numeric(3));
66     m3.set(1,0,numeric(3)).set(1,1,numeric(-2)).set(1,2,numeric(2));
67     m3.set(2,0,numeric(3)).set(2,1,numeric(4)).set(2,2,numeric(-2));
68     det = m3.determinant();
69     if (det != 42) {
70         clog << "determinant of 3x3 matrix " << m3
71              << " erroneously returned " << det << endl;
72         ++result;
73     }
74     
75     // check dense symbolic 2x2 matrix determinant
76     m2.set(0,0,a/(a-b)).set(0,1,numeric(1));
77     m2.set(1,0,b/(a-b)).set(1,1,numeric(1));
78     det = m2.determinant(true);
79     if (det != 1) {
80         clog << "determinant of 2x2 matrix " << m2
81              << " erroneously returned " << det << endl;
82         ++result;
83     }
84
85     // check characteristic polynomial
86     m3.set(0,0,a).set(0,1,-2).set(0,2,2);
87     m3.set(1,0,3).set(1,1,a-1).set(1,2,2);
88     m3.set(2,0,3).set(2,1,4).set(2,2,a-3);
89     ex p = m3.charpoly(a);
90     if (p != 0) {
91         clog << "charpoly of 3x3 matrix " << m3
92              << " erroneously returned " << p << endl;
93         ++result;
94     }
95     
96     return result;
97 }
98
99 static unsigned matrix_invert1(void)
100 {
101     matrix m(1,1);
102     symbol a("a");
103
104     m.set(0,0,a);
105     matrix m_i = m.inverse();
106     
107     if (m_i(0,0) != pow(a,-1)) {
108         clog << "inversion of 1x1 matrix " << m
109              << " erroneously returned " << m_i << endl;
110         return 1;
111     }
112     return 0;
113 }
114
115 static unsigned matrix_invert2(void)
116 {
117     matrix m(2,2);
118     symbol a("a"), b("b"), c("c"), d("d");
119     m.set(0,0,a).set(0,1,b);
120     m.set(1,0,c).set(1,1,d);
121     matrix m_i = m.inverse();
122     ex det = m.determinant().expand();
123     
124     if ( (normal(m_i(0,0)*det) != d) ||
125          (normal(m_i(0,1)*det) != -b) ||
126          (normal(m_i(1,0)*det) != -c) ||
127          (normal(m_i(1,1)*det) != a) ) {
128         clog << "inversion of 2x2 matrix " << m
129              << " erroneously returned " << m_i << endl;
130         return 1;
131     }
132     return 0;
133 }
134
135 static unsigned matrix_invert3(void)
136 {
137     matrix m(3,3);
138     symbol a("a"), b("b"), c("c");
139     symbol d("d"), e("e"), f("f");
140     symbol g("g"), h("h"), i("i");
141     m.set(0,0,a).set(0,1,b).set(0,2,c);
142     m.set(1,0,d).set(1,1,e).set(1,2,f);
143     m.set(2,0,g).set(2,1,h).set(2,2,i);
144     matrix m_i = m.inverse();
145     ex det = m.determinant().normal().expand();
146     
147     if ( (normal(m_i(0,0)*det) != (e*i-f*h)) ||
148          (normal(m_i(0,1)*det) != (c*h-b*i)) ||
149          (normal(m_i(0,2)*det) != (b*f-c*e)) ||
150          (normal(m_i(1,0)*det) != (f*g-d*i)) ||
151          (normal(m_i(1,1)*det) != (a*i-c*g)) ||
152          (normal(m_i(1,2)*det) != (c*d-a*f)) ||
153          (normal(m_i(2,0)*det) != (d*h-e*g)) ||
154          (normal(m_i(2,1)*det) != (b*g-a*h)) ||
155          (normal(m_i(2,2)*det) != (a*e-b*d)) ) {
156         clog << "inversion of 3x3 matrix " << m
157              << " erroneously returned " << m_i << endl;
158         return 1;
159     }
160     return 0;
161 }
162
163 static unsigned matrix_misc(void)
164 {
165     unsigned result = 0;
166     matrix m1(2,2);
167     symbol a("a"), b("b"), c("c"), d("d"), e("e"), f("f");
168     m1.set(0,0,a).set(0,1,b);
169     m1.set(1,0,c).set(1,1,d);
170     ex tr = trace(m1);
171     
172     // check a simple trace
173     if (tr.compare(a+d)) {
174         clog << "trace of 2x2 matrix " << m1
175              << " erroneously returned " << tr << endl;
176         ++result;
177     }
178     
179     // and two simple transpositions
180     matrix m2 = transpose(m1);
181     if (m2(0,0) != a || m2(0,1) != c || m2(1,0) != b || m2(1,1) != d) {
182         clog << "transpose of 2x2 matrix " << m1
183              << " erroneously returned " << m2 << endl;
184         ++result;
185     }
186     matrix m3(3,2);
187     m3.set(0,0,a).set(0,1,b);
188     m3.set(1,0,c).set(1,1,d);
189     m3.set(2,0,e).set(2,1,f);
190     if (transpose(transpose(m3)) != m3) {
191         clog << "transposing 3x2 matrix " << m3 << " twice"
192              << " erroneously returned " << transpose(transpose(m3)) << endl;
193         ++result;
194     }
195     
196     // produce a runtime-error by inverting a singular matrix and catch it
197     matrix m4(2,2);
198     matrix m5;
199     bool caught=false;
200     try {
201         m5 = inverse(m4);
202     }
203     catch (std::runtime_error err) {
204         caught=true;
205     }
206     if (!caught) {
207         cerr << "singular 2x2 matrix " << m4
208              << " erroneously inverted to " << m5 << endl;
209         ++result;
210     }
211     
212     return result;
213 }
214
215 unsigned matrix_checks(void)
216 {
217     unsigned result = 0;
218     
219     cout << "checking symbolic matrix manipulations..." << flush;
220     clog << "---------symbolic matrix manipulations:" << endl;
221     
222     result += matrix_determinants();
223     result += matrix_invert1();
224     result += matrix_invert2();
225     result += matrix_invert3();
226     result += matrix_misc();
227     
228     if (! result) {
229         cout << " passed ";
230         clog << "(no output)" << endl;
231     } else {
232         cout << " failed ";
233     }
234     
235     return result;
236 }