check/check_matrices.cpp

   1 /** @file check_matrices.cpp
   2  *
   3  *  Here we test manipulations on GiNaC's symbolic matrices.  They are a
   4  *  well-tried resource for cross-checking the underlying symbolic
   5  *  manipulations. */
   6
   7 /*
   8  *  GiNaC Copyright (C) 1999-2010 Johannes Gutenberg University Mainz, Germany
   9  *
  10  *  This program is free software; you can redistribute it and/or modify
  11  *  it under the terms of the GNU General Public License as published by
  12  *  the Free Software Foundation; either version 2 of the License, or
  13  *  (at your option) any later version.
  14  *
  15  *  This program is distributed in the hope that it will be useful,
  16  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  18  *  GNU General Public License for more details.
  19  *
  20  *  You should have received a copy of the GNU General Public License
  21  *  along with this program; if not, write to the Free Software
  22  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  23  */
  24
  25 #include "ginac.h"
  26 using namespace GiNaC;
  27
  28 #include <cstdlib> // for rand(), RAND_MAX
  29 #include <iostream>
  30 using namespace std;
  31
  32 extern const ex
  33 sparse_tree(const symbol & x, const symbol & y, const symbol & z,
  34             int level, bool trig = false, bool rational = true,
  35             bool complex = false);
  36 extern const ex
  37 dense_univariate_poly(const symbol & x, unsigned degree);
  38
  39 /* determinants of some sparse symbolic matrices with coefficients in
  40  * an integral domain. */
  41 static unsigned integdom_matrix_determinants()
  42 {
  43         unsigned result = 0;
  44         symbol a("a");
  45
  46         for (unsigned size=3; size<22; ++size) {
  47                 matrix A(size,size);
  48                 // populate one element in each row:
  49                 for (unsigned r=0; r<size-1; ++r)
  50                         A.set(r,unsigned(rand()%size),dense_univariate_poly(a,5));
  51                 // set the last row to a linear combination of two other lines
  52                 // to guarantee that the determinant is zero:
  53                 for (unsigned c=0; c<size; ++c)
  54                         A.set(size-1,c,A(0,c)-A(size-2,c));
  55                 if (!A.determinant().is_zero()) {
  56                         clog << "Determinant of " << size << "x" << size << " matrix "
  57                              << endl << A << endl
  58                              << "was not found to vanish!" << endl;
  59                         ++result;
  60                 }
  61         }
  62
  63         return result;
  64 }
  65
  66 /* determinants of some symbolic matrices with multivariate rational function
  67  * coefficients. */
  68 static unsigned rational_matrix_determinants()
  69 {
  70         unsigned result = 0;
  71         symbol a("a"), b("b"), c("c");
  72
  73         for (unsigned size=3; size<9; ++size) {
  74                 matrix A(size,size);
  75                 for (unsigned r=0; r<size-1; ++r) {
  76                         // populate one or two elements in each row:
  77                         for (unsigned ec=0; ec<2; ++ec) {
  78                                 ex numer = sparse_tree(a, b, c, 1+rand()%4, false, false, false);
  79                                 ex denom;
  80                                 do {
  81                                         denom = sparse_tree(a, b, c, rand()%2, false, false, false);
  82                                 } while (denom.is_zero());
  83                                 A.set(r,unsigned(rand()%size),numer/denom);
  84                         }
  85                 }
  86                 // set the last row to a linear combination of two other lines
  87                 // to guarantee that the determinant is zero:
  88                 for (unsigned co=0; co<size; ++co)
  89                         A.set(size-1,co,A(0,co)-A(size-2,co));
  90                 if (!A.determinant().is_zero()) {
  91                         clog << "Determinant of " << size << "x" << size << " matrix "
  92                              << endl << A << endl
  93                              << "was not found to vanish!" << endl;
  94                         ++result;
  95                 }
  96         }
  97
  98         return result;
  99 }
 100
 101 /* Some quite funny determinants with functions and stuff like that inside. */
 102 static unsigned funny_matrix_determinants()
 103 {
 104         unsigned result = 0;
 105         symbol a("a"), b("b"), c("c");
 106
 107         for (unsigned size=3; size<8; ++size) {
 108                 matrix A(size,size);
 109                 for (unsigned co=0; co<size-1; ++co) {
 110                         // populate one or two elements in each row:
 111                         for (unsigned ec=0; ec<2; ++ec) {
 112                                 ex numer = sparse_tree(a, b, c, 1+rand()%3, true, true, false);
 113                                 ex denom;
 114                                 do {
 115                                         denom = sparse_tree(a, b, c, rand()%2, false, true, false);
 116                                 } while (denom.is_zero());
 117                                 A.set(unsigned(rand()%size),co,numer/denom);
 118                         }
 119                 }
 120                 // set the last column to a linear combination of two other columns
 121                 // to guarantee that the determinant is zero:
 122                 for (unsigned ro=0; ro<size; ++ro)
 123                         A.set(ro,size-1,A(ro,0)-A(ro,size-2));
 124                 if (!A.determinant().is_zero()) {
 125                         clog << "Determinant of " << size << "x" << size << " matrix "
 126                              << endl << A << endl
 127                              << "was not found to vanish!" << endl;
 128                         ++result;
 129                 }
 130         }
 131
 132         return result;
 133 }
 134
 135 /* compare results from different determinant algorithms.*/
 136 static unsigned compare_matrix_determinants()
 137 {
 138         unsigned result = 0;
 139         symbol a("a");
 140
 141         for (unsigned size=2; size<8; ++size) {
 142                 matrix A(size,size);
 143                 for (unsigned co=0; co<size; ++co) {
 144                         for (unsigned ro=0; ro<size; ++ro) {
 145                                 // populate some elements
 146                                 ex elem = 0;
 147                                 if (rand()%(size/2) == 0)
 148                                         elem = sparse_tree(a, a, a, rand()%3, false, true, false);
 149                                 A.set(ro,co,elem);
 150                         }
 151                 }
 152                 ex det_gauss = A.determinant(determinant_algo::gauss);
 153                 ex det_laplace = A.determinant(determinant_algo::laplace);
 154                 ex det_divfree = A.determinant(determinant_algo::divfree);
 155                 ex det_bareiss = A.determinant(determinant_algo::bareiss);
 156                 if ((det_gauss-det_laplace).normal() != 0 ||
 157                         (det_bareiss-det_laplace).normal() != 0 ||
 158                         (det_divfree-det_laplace).normal() != 0) {
 159                         clog << "Determinant of " << size << "x" << size << " matrix "
 160                              << endl << A << endl
 161                              << "is inconsistent between different algorithms:" << endl
 162                              << "Gauss elimination:   " << det_gauss << endl
 163                              << "Minor elimination:   " << det_laplace << endl
 164                              << "Division-free elim.: " << det_divfree << endl
 165                              << "Fraction-free elim.: " << det_bareiss << endl;
 166                         ++result;
 167                 }
 168         }
 169
 170         return result;
 171 }
 172
 173 static unsigned symbolic_matrix_inverse()
 174 {
 175         unsigned result = 0;
 176         symbol a("a"), b("b"), c("c");
 177
 178         for (unsigned size=2; size<6; ++size) {
 179                 matrix A(size,size);
 180                 do {
 181                         for (unsigned co=0; co<size; ++co) {
 182                                 for (unsigned ro=0; ro<size; ++ro) {
 183                                         // populate some elements
 184                                         ex elem = 0;
 185                                         if (rand()%(size/2) == 0)
 186                                                 elem = sparse_tree(a, b, c, rand()%2, false, true, false);
 187                                         A.set(ro,co,elem);
 188                                 }
 189                         }
 190                 } while (A.determinant() == 0);
 191                 matrix B = A.inverse();
 192                 matrix C = A.mul(B);
 193                 bool ok = true;
 194                 for (unsigned ro=0; ro<size; ++ro)
 195                         for (unsigned co=0; co<size; ++co)
 196                                 if (C(ro,co).normal() != (ro==co?1:0))
 197                                         ok = false;
 198                 if (!ok) {
 199                         clog << "Inverse of " << size << "x" << size << " matrix "
 200                              << endl << A << endl
 201                              << "erroneously returned: "
 202                              << endl << B << endl;
 203                         ++result;
 204                 }
 205         }
 206
 207         return result;
 208 }
 209
 210 unsigned check_matrices()
 211 {
 212         unsigned result = 0;
 213
 214         cout << "checking symbolic matrix manipulations" << flush;
 215
 216         result += integdom_matrix_determinants();  cout << '.' << flush;
 217         result += rational_matrix_determinants();  cout << '.' << flush;
 218         result += funny_matrix_determinants();  cout << '.' << flush;
 219         result += compare_matrix_determinants();  cout << '.' << flush;
 220         result += symbolic_matrix_inverse();  cout << '.' << flush;
 221
 222         return result;
 223 }
 224
 225 int main(int argc, char** argv)
 226 {
 227         return check_matrices();
 228 }