check/matrix_checks.cpp

   1 /** @file matrix_checks.cpp
   2  *
   3  *  Here we test manipulations on GiNaC's symbolic matrices. */
   4
   5 /*
   6  *  GiNaC Copyright (C) 1999-2000 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  21  */
  22
  23 #include <stdexcept>
  24 #include "ginac.h"
  25
  26 #ifndef NO_NAMESPACE_GINAC
  27 using namespace GiNaC;
  28 #endif // ndef NO_NAMESPACE_GINAC
  29
  30 static unsigned matrix_determinants(void)
  31 {
  32     unsigned result = 0;
  33     ex det;
  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");
  38
  39     // check symbolic trivial matrix determinant
  40     m1.set(0,0,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.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;
  55         ++result;
  56     }
  57
  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;
  66         ++result;
  67     }
  68
  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();
  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.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);
  84     if (det != 1) {
  85         clog << "determinant of 2x2 matrix " << m2
  86              << " erroneously returned " << det << endl;
  87         ++result;
  88     }
  89
  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);
  95     if (p != 0) {
  96         clog << "charpoly of 3x3 matrix " << m3
  97              << " erroneously returned " << p << endl;
  98         ++result;
  99     }
 100
 101     return result;
 102 }
 103
 104 static unsigned matrix_invert1(void)
 105 {
 106     matrix m(1,1);
 107     symbol a("a");
 108
 109     m.set(0,0,a);
 110     matrix m_i = m.inverse();
 111
 112     if (m_i(0,0) != pow(a,-1)) {
 113         clog << "inversion of 1x1 matrix " << m
 114              << " erroneously returned " << m_i << endl;
 115         return 1;
 116     }
 117     return 0;
 118 }
 119
 120 static unsigned matrix_invert2(void)
 121 {
 122     matrix m(2,2);
 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();
 128
 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;
 135         return 1;
 136     }
 137     return 0;
 138 }
 139
 140 static unsigned matrix_invert3(void)
 141 {
 142     matrix m(3,3);
 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();
 151
 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;
 163         return 1;
 164     }
 165     return 0;
 166 }
 167
 168 static unsigned matrix_misc(void)
 169 {
 170     unsigned result = 0;
 171     matrix m1(2,2);
 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);
 175     ex tr = trace(m1);
 176
 177     // check a simple trace
 178     if (tr.compare(a+d)) {
 179         clog << "trace of 2x2 matrix " << m1
 180              << " erroneously returned " << tr << endl;
 181         ++result;
 182     }
 183
 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;
 189         ++result;
 190     }
 191     matrix m3(3,2);
 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;
 198         ++result;
 199     }
 200
 201     // produce a runtime-error by inverting a singular matrix and catch it
 202     matrix m4(2,2);
 203     matrix m5;
 204     bool caught=false;
 205     try {
 206         m5 = inverse(m4);
 207     }
 208     catch (std::runtime_error err) {
 209         caught=true;
 210     }
 211     if (!caught) {
 212         cerr << "singular 2x2 matrix " << m4
 213              << " erroneously inverted to " << m5 << endl;
 214         ++result;
 215     }
 216
 217     return result;
 218 }
 219
 220 unsigned matrix_checks(void)
 221 {
 222     unsigned result = 0;
 223
 224     cout << "checking symbolic matrix manipulations..." << flush;
 225     clog << "---------symbolic matrix manipulations:" << endl;
 226
 227     result += matrix_determinants();
 228     result += matrix_invert1();
 229     result += matrix_invert2();
 230     result += matrix_invert3();
 231     result += matrix_misc();
 232
 233     if (!result) {
 234         cout << " passed ";
 235         clog << "(no output)" << endl;
 236     } else {
 237         cout << " failed ";
 238     }
 239
 240     return result;
 241 }