1 // check/matrix_checks.cpp
3 /* Here we test manipulations on GiNaC's symbolic matrices. */
6 #include <GiNaC/ginac.h>
8 static unsigned matrix_determinants(void)
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
19 det = m1.determinant();
21 clog << "determinant of 1x1 matrix " << m1
22 << " erroneously returned " << det << endl;
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;
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;
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();
53 clog << "determinant of 3x3 matrix " << m3
54 << " erroneously returned " << det << endl;
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);
63 clog << "determinant of 2x2 matrix " << m2
64 << " erroneously returned " << det << endl;
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);
74 clog << "charpoly of 3x3 matrix " << m3
75 << " erroneously returned " << p << endl;
82 static unsigned matrix_invert1(void)
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;
98 static unsigned matrix_invert2(void)
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;
118 static unsigned matrix_invert3(void)
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;
146 static unsigned matrix_misc(void)
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);
155 // check a simple trace
156 if (tr.compare(a+d)) {
157 clog << "trace of 2x2 matrix " << m1
158 << " erroneously returned " << tr << endl;
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;
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;
179 // produce a runtime-error by inverting a singular matrix and catch it
186 catch (std::runtime_error err) {
190 cerr << "singular 2x2 matrix " << m4
191 << " erroneously inverted to " << m5 << endl;
198 unsigned matrix_checks(void)
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();
213 clog << "(no output)" << endl;