check/matrix_checks.cpp

   1 // check/matrix_checks.cpp
   2
   3 /* Here we test manipulations on GiNaC's symbolic matrices. */
   4
   5 #include <stdexcept>
   6 #include <ginac/ginac.h>
   7
   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");
  16
  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     }
  25
  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     }
  35
  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     }
  46
  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     }
  57
  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     }
  67
  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     }
  78
  79     return result;
  80 }
  81
  82 static unsigned matrix_invert1(void)
  83 {
  84     matrix m(1,1);
  85     symbol a("a");
  86
  87     m.set(0,0,a);
  88     matrix m_i = m.inverse();
  89
  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 }
  97
  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();
 106
 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 }
 117
 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();
 129
 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 }
 145
 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);
 154
 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     }
 161
 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     }
 178
 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     }
 194
 195     return result;
 196 }
 197
 198 unsigned matrix_checks(void)
 199 {
 200     unsigned result = 0;
 201
 202     cout << "checking symbolic matrix manipulations..." << flush;
 203     clog << "---------symbolic matrix manipulations:" << endl;
 204
 205     result += matrix_determinants();
 206     result += matrix_invert1();
 207     result += matrix_invert2();
 208     result += matrix_invert3();
 209     result += matrix_misc();
 210
 211     if (! result) {
 212         cout << " passed ";
 213         clog << "(no output)" << endl;
 214     } else {
 215         cout << " failed ";
 216     }
 217
 218     return result;
 219 }