586c52f57673a0bc20a33e89ba152f01697e31f7
[ginac.git] / check / exam_matrices.cpp
1 /** @file exam_matrices.cpp
2  *
3  *  Here we examine manipulations on GiNaC's symbolic matrices. */
4
5 /*
6  *  GiNaC Copyright (C) 1999-2018 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., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
21  */
22
23 #include "ginac.h"
24 using namespace GiNaC;
25
26 #include <iostream>
27 #include <stdexcept>
28 using namespace std;
29
30 static unsigned matrix_determinants()
31 {
32         unsigned result = 0;
33         ex det;
34         matrix m1(1,1), m2(2,2), m3(3,3), m4(4,4);
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 = matrix{{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 = matrix{{a, b},
50                     {c, 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 = matrix{{a, b, c},
60                     {d, e, f},
61                     {g, h, i}};
62         det = m3.determinant();
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 = matrix{{0, -1,  3},
71                     {3, -2,  2},
72                     {3,  4, -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 = matrix{{a/(a-b), 1},
82                     {b/(a-b), 1}};
83         det = m2.determinant();
84         if (det != 1) {
85                 if (det.normal() == 1)  // only half wrong
86                         clog << "determinant of 2x2 matrix " << m2
87                              << " was returned unnormalized as " << det << endl;
88                 else  // totally wrong
89                         clog << "determinant of 2x2 matrix " << m2
90                              << " erroneously returned " << det << endl;
91                 ++result;
92         }
93         
94         // check sparse symbolic 4x4 matrix determinant
95         m4.set(0,1,a).set(1,0,b).set(3,2,c).set(2,3,d);
96         det = m4.determinant();
97         if (det != a*b*c*d) {
98                 clog << "determinant of 4x4 matrix " << m4
99                      << " erroneously returned " << det << endl;
100                 ++result;
101         }
102         
103         // check characteristic polynomial
104         m3 = matrix{{a, -2,   2},
105                     {3, a-1,  2},
106                     {3,  4,  a-3}};
107         ex p = m3.charpoly(a);
108         if (p != 0) {
109                 clog << "charpoly of 3x3 matrix " << m3
110                      << " erroneously returned " << p << endl;
111                 ++result;
112         }
113         
114         return result;
115 }
116
117 static unsigned matrix_invert1()
118 {
119         unsigned result = 0;
120         matrix m(1,1);
121         symbol a("a");
122         
123         m.set(0,0,a);
124         matrix m_i = m.inverse();
125         
126         if (m_i(0,0) != pow(a,-1)) {
127                 clog << "inversion of 1x1 matrix " << m
128                      << " erroneously returned " << m_i << endl;
129                 ++result;
130         }
131         
132         return result;
133 }
134
135 static unsigned matrix_invert2()
136 {
137         unsigned result = 0;
138         symbol a("a"), b("b"), c("c"), d("d");
139         matrix m = {{a, b},
140                     {c, d}};
141         matrix m_i = m.inverse();
142         ex det = m.determinant();
143         
144         if ((normal(m_i(0,0)*det) != d) ||
145             (normal(m_i(0,1)*det) != -b) ||
146             (normal(m_i(1,0)*det) != -c) ||
147             (normal(m_i(1,1)*det) != a)) {
148                 clog << "inversion of 2x2 matrix " << m
149                      << " erroneously returned " << m_i << endl;
150                 ++result;
151         }
152         
153         return result;
154 }
155
156 static unsigned matrix_invert3()
157 {
158         unsigned result = 0;
159         symbol a("a"), b("b"), c("c");
160         symbol d("d"), e("e"), f("f");
161         symbol g("g"), h("h"), i("i");
162         matrix m = {{a, b, c},
163                     {d, e, f},
164                     {g, h, i}};
165         matrix m_i = m.inverse();
166         ex det = m.determinant();
167         
168         if ((normal(m_i(0,0)*det) != (e*i-f*h)) ||
169             (normal(m_i(0,1)*det) != (c*h-b*i)) ||
170             (normal(m_i(0,2)*det) != (b*f-c*e)) ||
171             (normal(m_i(1,0)*det) != (f*g-d*i)) ||
172             (normal(m_i(1,1)*det) != (a*i-c*g)) ||
173             (normal(m_i(1,2)*det) != (c*d-a*f)) ||
174             (normal(m_i(2,0)*det) != (d*h-e*g)) ||
175             (normal(m_i(2,1)*det) != (b*g-a*h)) ||
176             (normal(m_i(2,2)*det) != (a*e-b*d))) {
177                 clog << "inversion of 3x3 matrix " << m
178                      << " erroneously returned " << m_i << endl;
179                 ++result;
180         }
181         
182         return result;
183 }
184
185 static unsigned matrix_solve2()
186 {
187         // check the solution of the multiple system A*X = B:
188         //       [ 1  2 -1 ] [ x0 y0 ]   [ 4 0 ]
189         //       [ 1  4 -2 ]*[ x1 y1 ] = [ 7 0 ]
190         //       [ a -2  2 ] [ x2 y2 ]   [ a 4 ]
191         unsigned result = 0;
192         symbol a("a");
193         symbol x0("x0"), x1("x1"), x2("x2");
194         symbol y0("y0"), y1("y1"), y2("y2");
195         matrix A = {{1,  2, -1},
196                     {1,  4, -2},
197                     {a, -2,  2}};
198         matrix B = {{4, 0},
199                     {7, 0},
200                     {a, 4}};
201         matrix X = {{x0 ,y0},
202                     {x1, y1},
203                     {x2, y2}};
204         matrix cmp = {{1, 0},
205                       {3, 2},
206                       {3, 4}};
207         matrix sol(A.solve(X, B));
208         if (cmp != sol) {
209                 clog << "Solving " << A << " * " << X << " == " << B << endl
210                      << "erroneously returned " << sol << endl;
211                 result = 1;
212         }
213
214         return result;
215 }
216
217 static unsigned matrix_evalm()
218 {
219         unsigned result = 0;
220
221         matrix S {{1, 2},
222                   {3, 4}};
223         matrix T {{1, 1},
224                   {2, -1}};
225         matrix R {{27, 14},
226                   {36, 26}};
227
228         ex e = ((S + T) * (S + 2*T));
229         ex f = e.evalm();
230         if (!f.is_equal(R)) {
231                 clog << "Evaluating " << e << " erroneously returned " << f << " instead of " << R << endl;
232                 result++;
233         }
234
235         return result;
236 }
237
238 static unsigned matrix_rank()
239 {
240         unsigned result = 0;
241         symbol x("x"), y("y");
242         matrix m(3,3);
243
244         // the zero matrix always has rank 0
245         if (m.rank() != 0) {
246                 clog << "The rank of " << m << " was not computed correctly." << endl;
247                 ++result;
248         }
249
250         // a trivial rank one example
251         m = {{1, 0, 0},
252              {2, 0, 0},
253              {3, 0, 0}};
254         if (m.rank() != 1) {
255                 clog << "The rank of " << m << " was not computed correctly." << endl;
256                 ++result;
257         }
258
259         // an example from Maple's help with rank two
260         m = {{x,   1,  0},
261              {0,   0,  1},
262              {x*y, y,  1}};
263         if (m.rank() != 2) {
264                 clog << "The rank of " << m << " was not computed correctly." << endl;
265                 ++result;
266         }
267
268         // the 3x3 unit matrix has rank 3
269         m = ex_to<matrix>(unit_matrix(3,3));
270         if (m.rank() != 3) {
271                 clog << "The rank of " << m << " was not computed correctly." << endl;
272                 ++result;
273         }
274
275         return result;  
276 }
277
278 static unsigned matrix_misc()
279 {
280         unsigned result = 0;
281         symbol a("a"), b("b"), c("c"), d("d"), e("e"), f("f");
282         matrix m1 = {{a, b},
283                      {c, d}};
284         ex tr = trace(m1);
285         
286         // check a simple trace
287         if (tr.compare(a+d)) {
288                 clog << "trace of 2x2 matrix " << m1
289                      << " erroneously returned " << tr << endl;
290                 ++result;
291         }
292         
293         // and two simple transpositions
294         matrix m2 = transpose(m1);
295         if (m2(0,0) != a || m2(0,1) != c || m2(1,0) != b || m2(1,1) != d) {
296                 clog << "transpose of 2x2 matrix " << m1
297                          << " erroneously returned " << m2 << endl;
298                 ++result;
299         }
300         matrix m3 = {{a, b},
301                      {c, d},
302                      {e, f}};
303         if (transpose(transpose(m3)) != m3) {
304                 clog << "transposing 3x2 matrix " << m3 << " twice"
305                      << " erroneously returned " << transpose(transpose(m3)) << endl;
306                 ++result;
307         }
308         
309         // produce a runtime-error by inverting a singular matrix and catch it
310         matrix m4(2,2);
311         matrix m5;
312         bool caught = false;
313         try {
314                 m5 = inverse(m4);
315         } catch (std::runtime_error err) {
316                 caught = true;
317         }
318         if (!caught) {
319                 cerr << "singular 2x2 matrix " << m4
320                      << " erroneously inverted to " << m5 << endl;
321                 ++result;
322         }
323         
324         return result;
325 }
326
327 unsigned exam_matrices()
328 {
329         unsigned result = 0;
330         
331         cout << "examining symbolic matrix manipulations" << flush;
332         
333         result += matrix_determinants();  cout << '.' << flush;
334         result += matrix_invert1();  cout << '.' << flush;
335         result += matrix_invert2();  cout << '.' << flush;
336         result += matrix_invert3();  cout << '.' << flush;
337         result += matrix_solve2();  cout << '.' << flush;
338         result += matrix_evalm();  cout << "." << flush;
339         result += matrix_rank();  cout << "." << flush;
340         result += matrix_misc();  cout << '.' << flush;
341         
342         return result;
343 }
344
345 int main(int argc, char** argv)
346 {
347         return exam_matrices();
348 }