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 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/ginac.h>
  25 using namespace GiNaC;
  26
  27 static unsigned matrix_determinants(void)
  28 {
  29     unsigned result = 0;
  30     ex det;
  31     matrix m1(1,1), m2(2,2), m3(3,3);
  32     symbol a("a"), b("b"), c("c");
  33     symbol d("d"), e("e"), f("f");
  34     symbol g("g"), h("h"), i("i");
  35
  36     // check symbolic trivial matrix determinant
  37     m1.set(0,0,a);
  38     det = m1.determinant();
  39     if (det != a) {
  40         clog << "determinant of 1x1 matrix " << m1
  41              << " erroneously returned " << det << endl;
  42         ++result;
  43     }
  44
  45     // check generic dense symbolic 2x2 matrix determinant
  46     m2.set(0,0,a).set(0,1,b);
  47     m2.set(1,0,c).set(1,1,d);
  48     det = m2.determinant();
  49     if (det != (a*d-b*c)) {
  50         clog << "determinant of 2x2 matrix " << m2
  51              << " erroneously returned " << det << endl;
  52         ++result;
  53     }
  54
  55     // check generic dense symbolic 3x3 matrix determinant
  56     m3.set(0,0,a).set(0,1,b).set(0,2,c);
  57     m3.set(1,0,d).set(1,1,e).set(1,2,f);
  58     m3.set(2,0,g).set(2,1,h).set(2,2,i);
  59     det = m3.determinant().expand();
  60     if (det != (a*e*i - a*f*h - d*b*i + d*c*h + g*b*f - g*c*e)) {
  61         clog << "determinant of 3x3 matrix " << m3
  62              << " erroneously returned " << det << endl;
  63         ++result;
  64     }
  65
  66     // check dense numeric 3x3 matrix determinant
  67     m3.set(0,0,numeric(0)).set(0,1,numeric(-1)).set(0,2,numeric(3));
  68     m3.set(1,0,numeric(3)).set(1,1,numeric(-2)).set(1,2,numeric(2));
  69     m3.set(2,0,numeric(3)).set(2,1,numeric(4)).set(2,2,numeric(-2));
  70     det = m3.determinant();
  71     if (det != 42) {
  72         clog << "determinant of 3x3 matrix " << m3
  73              << " erroneously returned " << det << endl;
  74         ++result;
  75     }
  76
  77     // check dense symbolic 2x2 matrix determinant
  78     m2.set(0,0,a/(a-b)).set(0,1,numeric(1));
  79     m2.set(1,0,b/(a-b)).set(1,1,numeric(1));
  80     det = m2.determinant(true);
  81     if (det != 1) {
  82         clog << "determinant of 2x2 matrix " << m2
  83              << " erroneously returned " << det << endl;
  84         ++result;
  85     }
  86
  87     // check characteristic polynomial
  88     m3.set(0,0,a).set(0,1,-2).set(0,2,2);
  89     m3.set(1,0,3).set(1,1,a-1).set(1,2,2);
  90     m3.set(2,0,3).set(2,1,4).set(2,2,a-3);
  91     ex p = m3.charpoly(a);
  92     if (p != 0) {
  93         clog << "charpoly of 3x3 matrix " << m3
  94              << " erroneously returned " << p << endl;
  95         ++result;
  96     }
  97
  98     return result;
  99 }
 100
 101 static unsigned matrix_invert1(void)
 102 {
 103     matrix m(1,1);
 104     symbol a("a");
 105
 106     m.set(0,0,a);
 107     matrix m_i = m.inverse();
 108
 109     if (m_i(0,0) != pow(a,-1)) {
 110         clog << "inversion of 1x1 matrix " << m
 111              << " erroneously returned " << m_i << endl;
 112         return 1;
 113     }
 114     return 0;
 115 }
 116
 117 static unsigned matrix_invert2(void)
 118 {
 119     matrix m(2,2);
 120     symbol a("a"), b("b"), c("c"), d("d");
 121     m.set(0,0,a).set(0,1,b);
 122     m.set(1,0,c).set(1,1,d);
 123     matrix m_i = m.inverse();
 124     ex det = m.determinant().expand();
 125
 126     if ( (normal(m_i(0,0)*det) != d) ||
 127          (normal(m_i(0,1)*det) != -b) ||
 128          (normal(m_i(1,0)*det) != -c) ||
 129          (normal(m_i(1,1)*det) != a) ) {
 130         clog << "inversion of 2x2 matrix " << m
 131              << " erroneously returned " << m_i << endl;
 132         return 1;
 133     }
 134     return 0;
 135 }
 136
 137 static unsigned matrix_invert3(void)
 138 {
 139     matrix m(3,3);
 140     symbol a("a"), b("b"), c("c");
 141     symbol d("d"), e("e"), f("f");
 142     symbol g("g"), h("h"), i("i");
 143     m.set(0,0,a).set(0,1,b).set(0,2,c);
 144     m.set(1,0,d).set(1,1,e).set(1,2,f);
 145     m.set(2,0,g).set(2,1,h).set(2,2,i);
 146     matrix m_i = m.inverse();
 147     ex det = m.determinant().normal().expand();
 148
 149     if ( (normal(m_i(0,0)*det) != (e*i-f*h)) ||
 150          (normal(m_i(0,1)*det) != (c*h-b*i)) ||
 151          (normal(m_i(0,2)*det) != (b*f-c*e)) ||
 152          (normal(m_i(1,0)*det) != (f*g-d*i)) ||
 153          (normal(m_i(1,1)*det) != (a*i-c*g)) ||
 154          (normal(m_i(1,2)*det) != (c*d-a*f)) ||
 155          (normal(m_i(2,0)*det) != (d*h-e*g)) ||
 156          (normal(m_i(2,1)*det) != (b*g-a*h)) ||
 157          (normal(m_i(2,2)*det) != (a*e-b*d)) ) {
 158         clog << "inversion of 3x3 matrix " << m
 159              << " erroneously returned " << m_i << endl;
 160         return 1;
 161     }
 162     return 0;
 163 }
 164
 165 static unsigned matrix_misc(void)
 166 {
 167     unsigned result = 0;
 168     matrix m1(2,2);
 169     symbol a("a"), b("b"), c("c"), d("d"), e("e"), f("f");
 170     m1.set(0,0,a).set(0,1,b);
 171     m1.set(1,0,c).set(1,1,d);
 172     ex tr = trace(m1);
 173
 174     // check a simple trace
 175     if (tr.compare(a+d)) {
 176         clog << "trace of 2x2 matrix " << m1
 177              << " erroneously returned " << tr << endl;
 178         ++result;
 179     }
 180
 181     // and two simple transpositions
 182     matrix m2 = transpose(m1);
 183     if (m2(0,0) != a || m2(0,1) != c || m2(1,0) != b || m2(1,1) != d) {
 184         clog << "transpose of 2x2 matrix " << m1
 185              << " erroneously returned " << m2 << endl;
 186         ++result;
 187     }
 188     matrix m3(3,2);
 189     m3.set(0,0,a).set(0,1,b);
 190     m3.set(1,0,c).set(1,1,d);
 191     m3.set(2,0,e).set(2,1,f);
 192     if (transpose(transpose(m3)) != m3) {
 193         clog << "transposing 3x2 matrix " << m3 << " twice"
 194              << " erroneously returned " << transpose(transpose(m3)) << endl;
 195         ++result;
 196     }
 197
 198     // produce a runtime-error by inverting a singular matrix and catch it
 199     matrix m4(2,2);
 200     matrix m5;
 201     bool caught=false;
 202     try {
 203         m5 = inverse(m4);
 204     }
 205     catch (std::runtime_error err) {
 206         caught=true;
 207     }
 208     if (!caught) {
 209         cerr << "singular 2x2 matrix " << m4
 210              << " erroneously inverted to " << m5 << endl;
 211         ++result;
 212     }
 213
 214     return result;
 215 }
 216
 217 unsigned matrix_checks(void)
 218 {
 219     unsigned result = 0;
 220
 221     cout << "checking symbolic matrix manipulations..." << flush;
 222     clog << "---------symbolic matrix manipulations:" << endl;
 223
 224     result += matrix_determinants();
 225     result += matrix_invert1();
 226     result += matrix_invert2();
 227     result += matrix_invert3();
 228     result += matrix_misc();
 229
 230     if (! result) {
 231         cout << " passed ";
 232         clog << "(no output)" << endl;
 233     } else {
 234         cout << " failed ";
 235     }
 236
 237     return result;
 238 }