1 // check/matrix_checks.cpp
3 /* Here we test manipulations on GiNaC's symbolic matrices. */
5 #include "ginac.h"
6 #include <stdexcept>
8 static unsigned matrix_determinants(void)
9 {
10     unsigned result = 0;
11     ex det;
12     matrix m1(1,1), m2(2,2), m3(3,3);
13     symbol a("a"), b("b"), c("c");
14     symbol d("d"), e("e"), f("f");
15     symbol g("g"), h("h"), i("i");
17     // check symbolic trivial matrix determinant
18     m1.set(0,0,a);
19     det = m1.determinant();
20     if (det != a) {
21         clog << "determinant of 1x1 matrix " << m1
22              << " erroneously returned " << det << endl;
23         ++result;
24     }
26     // check generic dense symbolic 2x2 matrix determinant
27     m2.set(0,0,a).set(0,1,b);
28     m2.set(1,0,c).set(1,1,d);
29     det = m2.determinant();
30     if (det != (a*d-b*c)) {
31         clog << "determinant of 2x2 matrix " << m2
32              << " erroneously returned " << det << endl;
33         ++result;
34     }
36     // check generic dense symbolic 3x3 matrix determinant
37     m3.set(0,0,a).set(0,1,b).set(0,2,c);
38     m3.set(1,0,d).set(1,1,e).set(1,2,f);
39     m3.set(2,0,g).set(2,1,h).set(2,2,i);
40     det = m3.determinant().expand();
41     if (det != (a*e*i - a*f*h - d*b*i + d*c*h + g*b*f - g*c*e)) {
42         clog << "determinant of 3x3 matrix " << m3
43              << " erroneously returned " << det << endl;
44         ++result;
45     }
47     // check dense numeric 3x3 matrix determinant
48     m3.set(0,0,numeric(0)).set(0,1,numeric(-1)).set(0,2,numeric(3));
49     m3.set(1,0,numeric(3)).set(1,1,numeric(-2)).set(1,2,numeric(2));
50     m3.set(2,0,numeric(3)).set(2,1,numeric(4)).set(2,2,numeric(-2));
51     det = m3.determinant();
52     if (det != 42) {
53         clog << "determinant of 3x3 matrix " << m3
54              << " erroneously returned " << det << endl;
55         ++result;
56     }
58     // check dense symbolic 2x2 matrix determinant
59     m2.set(0,0,a/(a-b)).set(0,1,numeric(1));
60     m2.set(1,0,b/(a-b)).set(1,1,numeric(1));
61     det = m2.determinant(true);
62     if (det != 1) {
63         clog << "determinant of 2x2 matrix " << m2
64              << " erroneously returned " << det << endl;
65         ++result;
66     }
68     // check characteristic polynomial
69     m3.set(0,0,a).set(0,1,-2).set(0,2,2);
70     m3.set(1,0,3).set(1,1,a-1).set(1,2,2);
71     m3.set(2,0,3).set(2,1,4).set(2,2,a-3);
72     ex p = m3.charpoly(a);
73     if (p != 0) {
74         clog << "charpoly of 3x3 matrix " << m3
75              << " erroneously returned " << p << endl;
76         ++result;
77     }
79     return result;
80 }
82 static unsigned matrix_invert1(void)
83 {
84     matrix m(1,1);
85     symbol a("a");
87     m.set(0,0,a);
88     matrix m_i = m.inverse();
90     if (m_i(0,0) != pow(a,-1)) {
91         clog << "inversion of 1x1 matrix " << m
92              << " erroneously returned " << m_i << endl;
93         return 1;
94     }
95     return 0;
96 }
98 static unsigned matrix_invert2(void)
99 {
100     matrix m(2,2);
101     symbol a("a"), b("b"), c("c"), d("d");
102     m.set(0,0,a).set(0,1,b);
103     m.set(1,0,c).set(1,1,d);
104     matrix m_i = m.inverse();
105     ex det = m.determinant().expand();
107     if ( (normal(m_i(0,0)*det) != d) ||
108          (normal(m_i(0,1)*det) != -b) ||
109          (normal(m_i(1,0)*det) != -c) ||
110          (normal(m_i(1,1)*det) != a) ) {
111         clog << "inversion of 2x2 matrix " << m
112              << " erroneously returned " << m_i << endl;
113         return 1;
114     }
115     return 0;
116 }
118 static unsigned matrix_invert3(void)
119 {
120     matrix m(3,3);
121     symbol a("a"), b("b"), c("c");
122     symbol d("d"), e("e"), f("f");
123     symbol g("g"), h("h"), i("i");
124     m.set(0,0,a).set(0,1,b).set(0,2,c);
125     m.set(1,0,d).set(1,1,e).set(1,2,f);
126     m.set(2,0,g).set(2,1,h).set(2,2,i);
127     matrix m_i = m.inverse();
128     ex det = m.determinant().normal().expand();
130     if ( (normal(m_i(0,0)*det) != (e*i-f*h)) ||
131          (normal(m_i(0,1)*det) != (c*h-b*i)) ||
132          (normal(m_i(0,2)*det) != (b*f-c*e)) ||
133          (normal(m_i(1,0)*det) != (f*g-d*i)) ||
134          (normal(m_i(1,1)*det) != (a*i-c*g)) ||
135          (normal(m_i(1,2)*det) != (c*d-a*f)) ||
136          (normal(m_i(2,0)*det) != (d*h-e*g)) ||
137          (normal(m_i(2,1)*det) != (b*g-a*h)) ||
138          (normal(m_i(2,2)*det) != (a*e-b*d)) ) {
139         clog << "inversion of 3x3 matrix " << m
140              << " erroneously returned " << m_i << endl;
141         return 1;
142     }
143     return 0;
144 }
146 static unsigned matrix_misc(void)
147 {
148     unsigned result = 0;
149     matrix m1(2,2);
150     symbol a("a"), b("b"), c("c"), d("d"), e("e"), f("f");
151     m1.set(0,0,a).set(0,1,b);
152     m1.set(1,0,c).set(1,1,d);
153     ex tr = trace(m1);
155     // check a simple trace
156     if (tr.compare(a+d)) {
157         clog << "trace of 2x2 matrix " << m1
158              << " erroneously returned " << tr << endl;
159         ++result;
160     }
162     // and two simple transpositions
163     matrix m2 = transpose(m1);
164     if (m2(0,0) != a || m2(0,1) != c || m2(1,0) != b || m2(1,1) != d) {
165         clog << "transpose of 2x2 matrix " << m1
166              << " erroneously returned " << m2 << endl;
167         ++result;
168     }
169     matrix m3(3,2);
170     m3.set(0,0,a).set(0,1,b);
171     m3.set(1,0,c).set(1,1,d);
172     m3.set(2,0,e).set(2,1,f);
173     if (transpose(transpose(m3)) != m3) {
174         clog << "transposing 3x2 matrix " << m3 << " twice"
175              << " erroneously returned " << transpose(transpose(m3)) << endl;
176         ++result;
177     }
179     // produce a runtime-error by inverting a singular matrix and catch it
180     matrix m4(2,2);
181     matrix m5;
182     bool caught=false;
183     try {
184         m5 = inverse(m4);
185     }
186     catch (std::runtime_error err) {
187         caught=true;
188     }
189     if (!caught) {
190         cerr << "singular 2x2 matrix " << m4
191              << " erroneously inverted to " << m5 << endl;
192         ++result;
193     }
195     return result;
196 }
198 unsigned matrix_checks(void)
199 {
200     unsigned result = 0;
202     cout << "checking symbolic matrix manipulations..." << flush;
203     clog << "---------symbolic matrix manipulations:" << endl;
205     result += matrix_determinants();
206     result += matrix_invert1();
207     result += matrix_invert2();
208     result += matrix_invert3();
209     result += matrix_misc();
211     if (! result) {
212         cout << " passed ";
213         clog << "(no output)" << endl;
214     } else {
215         cout << " failed ";
216     }
218     return result;
219 }