1 /** @file matrix_checks.cpp
3 * Here we test manipulations on GiNaC's symbolic matrices. */
6 * GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany
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.
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.
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
24 #include <ginac/ginac.h>
26 #ifndef NO_NAMESPACE_GINAC
27 using namespace GiNaC;
28 #endif // ndef NO_NAMESPACE_GINAC
30 static unsigned matrix_determinants(void)
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");
39 // check symbolic trivial matrix determinant
41 det = m1.determinant();
43 clog << "determinant of 1x1 matrix " << m1
44 << " erroneously returned " << det << endl;
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;
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;
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();
75 clog << "determinant of 3x3 matrix " << m3
76 << " erroneously returned " << det << endl;
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);
85 clog << "determinant of 2x2 matrix " << m2
86 << " erroneously returned " << det << endl;
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);
96 clog << "charpoly of 3x3 matrix " << m3
97 << " erroneously returned " << p << endl;
104 static unsigned matrix_invert1(void)
110 matrix m_i = m.inverse();
112 if (m_i(0,0) != pow(a,-1)) {
113 clog << "inversion of 1x1 matrix " << m
114 << " erroneously returned " << m_i << endl;
120 static unsigned matrix_invert2(void)
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();
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;
140 static unsigned matrix_invert3(void)
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();
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;
168 static unsigned matrix_misc(void)
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);
177 // check a simple trace
178 if (tr.compare(a+d)) {
179 clog << "trace of 2x2 matrix " << m1
180 << " erroneously returned " << tr << endl;
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;
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;
201 // produce a runtime-error by inverting a singular matrix and catch it
208 catch (std::runtime_error err) {
212 cerr << "singular 2x2 matrix " << m4
213 << " erroneously inverted to " << m5 << endl;
220 unsigned matrix_checks(void)
224 cout << "checking symbolic matrix manipulations..." << flush;
225 clog << "---------symbolic matrix manipulations:" << endl;
227 result += matrix_determinants();
228 result += matrix_invert1();
229 result += matrix_invert2();
230 result += matrix_invert3();
231 result += matrix_misc();
235 clog << "(no output)" << endl;