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-2015 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 }