check/matrix_checks.cpp

   1 /** @file matrix_checks.cpp
   2  *
   3  *  Here we test manipulations on GiNaC's symbolic matrices.
   4  *
   5  *  GiNaC Copyright (C) 1999 Johannes Gutenberg University Mainz, Germany
   6  *
   7  *  This program is free software; you can redistribute it and/or modify
   8  *  it under the terms of the GNU General Public License as published by
   9  *  the Free Software Foundation; either version 2 of the License, or
  10  *  (at your option) any later version.
  11  *
  12  *  This program is distributed in the hope that it will be useful,
  13  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15  *  GNU General Public License for more details.
  16  *
  17  *  You should have received a copy of the GNU General Public License
  18  *  along with this program; if not, write to the Free Software
  19  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  20  */
  21
  22 #include <stdexcept>
  23 #include <ginac/ginac.h>
  24
  25 static unsigned matrix_determinants(void)
  26 {
  27     unsigned result = 0;
  28     ex det;
  29     matrix m1(1,1), m2(2,2), m3(3,3);
  30     symbol a("a"), b("b"), c("c");
  31     symbol d("d"), e("e"), f("f");
  32     symbol g("g"), h("h"), i("i");
  33
  34     // check symbolic trivial matrix determinant
  35     m1.set(0,0,a);
  36     det = m1.determinant();
  37     if (det != a) {
  38         clog << "determinant of 1x1 matrix " << m1
  39              << " erroneously returned " << det << endl;
  40         ++result;
  41     }
  42
  43     // check generic dense symbolic 2x2 matrix determinant
  44     m2.set(0,0,a).set(0,1,b);
  45     m2.set(1,0,c).set(1,1,d);
  46     det = m2.determinant();
  47     if (det != (a*d-b*c)) {
  48         clog << "determinant of 2x2 matrix " << m2
  49              << " erroneously returned " << det << endl;
  50         ++result;
  51     }
  52
  53     // check generic dense symbolic 3x3 matrix determinant
  54     m3.set(0,0,a).set(0,1,b).set(0,2,c);
  55     m3.set(1,0,d).set(1,1,e).set(1,2,f);
  56     m3.set(2,0,g).set(2,1,h).set(2,2,i);
  57     det = m3.determinant().expand();
  58     if (det != (a*e*i - a*f*h - d*b*i + d*c*h + g*b*f - g*c*e)) {
  59         clog << "determinant of 3x3 matrix " << m3
  60              << " erroneously returned " << det << endl;
  61         ++result;
  62     }
  63
  64     // check dense numeric 3x3 matrix determinant
  65     m3.set(0,0,numeric(0)).set(0,1,numeric(-1)).set(0,2,numeric(3));
  66     m3.set(1,0,numeric(3)).set(1,1,numeric(-2)).set(1,2,numeric(2));
  67     m3.set(2,0,numeric(3)).set(2,1,numeric(4)).set(2,2,numeric(-2));
  68     det = m3.determinant();
  69     if (det != 42) {
  70         clog << "determinant of 3x3 matrix " << m3
  71              << " erroneously returned " << det << endl;
  72         ++result;
  73     }
  74
  75     // check dense symbolic 2x2 matrix determinant
  76     m2.set(0,0,a/(a-b)).set(0,1,numeric(1));
  77     m2.set(1,0,b/(a-b)).set(1,1,numeric(1));
  78     det = m2.determinant(true);
  79     if (det != 1) {
  80         clog << "determinant of 2x2 matrix " << m2
  81              << " erroneously returned " << det << endl;
  82         ++result;
  83     }
  84
  85     // check characteristic polynomial
  86     m3.set(0,0,a).set(0,1,-2).set(0,2,2);
  87     m3.set(1,0,3).set(1,1,a-1).set(1,2,2);
  88     m3.set(2,0,3).set(2,1,4).set(2,2,a-3);
  89     ex p = m3.charpoly(a);
  90     if (p != 0) {
  91         clog << "charpoly of 3x3 matrix " << m3
  92              << " erroneously returned " << p << endl;
  93         ++result;
  94     }
  95
  96     return result;
  97 }
  98
  99 static unsigned matrix_invert1(void)
 100 {
 101     matrix m(1,1);
 102     symbol a("a");
 103
 104     m.set(0,0,a);
 105     matrix m_i = m.inverse();
 106
 107     if (m_i(0,0) != pow(a,-1)) {
 108         clog << "inversion of 1x1 matrix " << m
 109              << " erroneously returned " << m_i << endl;
 110         return 1;
 111     }
 112     return 0;
 113 }
 114
 115 static unsigned matrix_invert2(void)
 116 {
 117     matrix m(2,2);
 118     symbol a("a"), b("b"), c("c"), d("d");
 119     m.set(0,0,a).set(0,1,b);
 120     m.set(1,0,c).set(1,1,d);
 121     matrix m_i = m.inverse();
 122     ex det = m.determinant().expand();
 123
 124     if ( (normal(m_i(0,0)*det) != d) ||
 125          (normal(m_i(0,1)*det) != -b) ||
 126          (normal(m_i(1,0)*det) != -c) ||
 127          (normal(m_i(1,1)*det) != a) ) {
 128         clog << "inversion of 2x2 matrix " << m
 129              << " erroneously returned " << m_i << endl;
 130         return 1;
 131     }
 132     return 0;
 133 }
 134
 135 static unsigned matrix_invert3(void)
 136 {
 137     matrix m(3,3);
 138     symbol a("a"), b("b"), c("c");
 139     symbol d("d"), e("e"), f("f");
 140     symbol g("g"), h("h"), i("i");
 141     m.set(0,0,a).set(0,1,b).set(0,2,c);
 142     m.set(1,0,d).set(1,1,e).set(1,2,f);
 143     m.set(2,0,g).set(2,1,h).set(2,2,i);
 144     matrix m_i = m.inverse();
 145     ex det = m.determinant().normal().expand();
 146
 147     if ( (normal(m_i(0,0)*det) != (e*i-f*h)) ||
 148          (normal(m_i(0,1)*det) != (c*h-b*i)) ||
 149          (normal(m_i(0,2)*det) != (b*f-c*e)) ||
 150          (normal(m_i(1,0)*det) != (f*g-d*i)) ||
 151          (normal(m_i(1,1)*det) != (a*i-c*g)) ||
 152          (normal(m_i(1,2)*det) != (c*d-a*f)) ||
 153          (normal(m_i(2,0)*det) != (d*h-e*g)) ||
 154          (normal(m_i(2,1)*det) != (b*g-a*h)) ||
 155          (normal(m_i(2,2)*det) != (a*e-b*d)) ) {
 156         clog << "inversion of 3x3 matrix " << m
 157              << " erroneously returned " << m_i << endl;
 158         return 1;
 159     }
 160     return 0;
 161 }
 162
 163 static unsigned matrix_misc(void)
 164 {
 165     unsigned result = 0;
 166     matrix m1(2,2);
 167     symbol a("a"), b("b"), c("c"), d("d"), e("e"), f("f");
 168     m1.set(0,0,a).set(0,1,b);
 169     m1.set(1,0,c).set(1,1,d);
 170     ex tr = trace(m1);
 171
 172     // check a simple trace
 173     if (tr.compare(a+d)) {
 174         clog << "trace of 2x2 matrix " << m1
 175              << " erroneously returned " << tr << endl;
 176         ++result;
 177     }
 178
 179     // and two simple transpositions
 180     matrix m2 = transpose(m1);
 181     if (m2(0,0) != a || m2(0,1) != c || m2(1,0) != b || m2(1,1) != d) {
 182         clog << "transpose of 2x2 matrix " << m1
 183              << " erroneously returned " << m2 << endl;
 184         ++result;
 185     }
 186     matrix m3(3,2);
 187     m3.set(0,0,a).set(0,1,b);
 188     m3.set(1,0,c).set(1,1,d);
 189     m3.set(2,0,e).set(2,1,f);
 190     if (transpose(transpose(m3)) != m3) {
 191         clog << "transposing 3x2 matrix " << m3 << " twice"
 192              << " erroneously returned " << transpose(transpose(m3)) << endl;
 193         ++result;
 194     }
 195
 196     // produce a runtime-error by inverting a singular matrix and catch it
 197     matrix m4(2,2);
 198     matrix m5;
 199     bool caught=false;
 200     try {
 201         m5 = inverse(m4);
 202     }
 203     catch (std::runtime_error err) {
 204         caught=true;
 205     }
 206     if (!caught) {
 207         cerr << "singular 2x2 matrix " << m4
 208              << " erroneously inverted to " << m5 << endl;
 209         ++result;
 210     }
 211
 212     return result;
 213 }
 214
 215 unsigned matrix_checks(void)
 216 {
 217     unsigned result = 0;
 218
 219     cout << "checking symbolic matrix manipulations..." << flush;
 220     clog << "---------symbolic matrix manipulations:" << endl;
 221
 222     result += matrix_determinants();
 223     result += matrix_invert1();
 224     result += matrix_invert2();
 225     result += matrix_invert3();
 226     result += matrix_misc();
 227
 228     if (! result) {
 229         cout << " passed ";
 230         clog << "(no output)" << endl;
 231     } else {
 232         cout << " failed ";
 233     }
 234
 235     return result;
 236 }