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-2008 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 <iostream>
  26 #include <cstdlib> // rand(), RAND_MAX
  27 #include "ginac.h"
  28 using namespace std;
  29 using namespace GiNaC;
  30
  31 extern const ex
  32 sparse_tree(const symbol & x, const symbol & y, const symbol & z,
  33             int level, bool trig = false, bool rational = true,
  34             bool complex = false);
  35 extern const ex
  36 dense_univariate_poly(const symbol & x, unsigned degree);
  37
  38 /* determinants of some sparse symbolic matrices with coefficients in
  39  * an integral domain. */
  40 static unsigned integdom_matrix_determinants()
  41 {
  42         unsigned result = 0;
  43         symbol a("a");
  44
  45         for (unsigned size=3; size<22; ++size) {
  46                 matrix A(size,size);
  47                 // populate one element in each row:
  48                 for (unsigned r=0; r<size-1; ++r)
  49                         A.set(r,unsigned(rand()%size),dense_univariate_poly(a,5));
  50                 // set the last row to a linear combination of two other lines
  51                 // to guarantee that the determinant is zero:
  52                 for (unsigned c=0; c<size; ++c)
  53                         A.set(size-1,c,A(0,c)-A(size-2,c));
  54                 if (!A.determinant().is_zero()) {
  55                         clog << "Determinant of " << size << "x" << size << " matrix "
  56                              << endl << A << endl
  57                              << "was not found to vanish!" << endl;
  58                         ++result;
  59                 }
  60         }
  61
  62         return result;
  63 }
  64
  65 /* determinants of some symbolic matrices with multivariate rational function
  66  * coefficients. */
  67 static unsigned rational_matrix_determinants()
  68 {
  69         unsigned result = 0;
  70         symbol a("a"), b("b"), c("c");
  71
  72         for (unsigned size=3; size<9; ++size) {
  73                 matrix A(size,size);
  74                 for (unsigned r=0; r<size-1; ++r) {
  75                         // populate one or two elements in each row:
  76                         for (unsigned ec=0; ec<2; ++ec) {
  77                                 ex numer = sparse_tree(a, b, c, 1+rand()%4, false, false, false);
  78                                 ex denom;
  79                                 do {
  80                                         denom = sparse_tree(a, b, c, rand()%2, false, false, false);
  81                                 } while (denom.is_zero());
  82                                 A.set(r,unsigned(rand()%size),numer/denom);
  83                         }
  84                 }
  85                 // set the last row to a linear combination of two other lines
  86                 // to guarantee that the determinant is zero:
  87                 for (unsigned co=0; co<size; ++co)
  88                         A.set(size-1,co,A(0,co)-A(size-2,co));
  89                 if (!A.determinant().is_zero()) {
  90                         clog << "Determinant of " << size << "x" << size << " matrix "
  91                              << endl << A << endl
  92                              << "was not found to vanish!" << endl;
  93                         ++result;
  94                 }
  95         }
  96
  97         return result;
  98 }
  99
 100 /* Some quite funny determinants with functions and stuff like that inside. */
 101 static unsigned funny_matrix_determinants()
 102 {
 103         unsigned result = 0;
 104         symbol a("a"), b("b"), c("c");
 105
 106         for (unsigned size=3; size<8; ++size) {
 107                 matrix A(size,size);
 108                 for (unsigned co=0; co<size-1; ++co) {
 109                         // populate one or two elements in each row:
 110                         for (unsigned ec=0; ec<2; ++ec) {
 111                                 ex numer = sparse_tree(a, b, c, 1+rand()%3, true, true, false);
 112                                 ex denom;
 113                                 do {
 114                                         denom = sparse_tree(a, b, c, rand()%2, false, true, false);
 115                                 } while (denom.is_zero());
 116                                 A.set(unsigned(rand()%size),co,numer/denom);
 117                         }
 118                 }
 119                 // set the last column to a linear combination of two other columns
 120                 // to guarantee that the determinant is zero:
 121                 for (unsigned ro=0; ro<size; ++ro)
 122                         A.set(ro,size-1,A(ro,0)-A(ro,size-2));
 123                 if (!A.determinant().is_zero()) {
 124                         clog << "Determinant of " << size << "x" << size << " matrix "
 125                              << endl << A << endl
 126                              << "was not found to vanish!" << endl;
 127                         ++result;
 128                 }
 129         }
 130
 131         return result;
 132 }
 133
 134 /* compare results from different determinant algorithms.*/
 135 static unsigned compare_matrix_determinants()
 136 {
 137         unsigned result = 0;
 138         symbol a("a");
 139
 140         for (unsigned size=2; size<8; ++size) {
 141                 matrix A(size,size);
 142                 for (unsigned co=0; co<size; ++co) {
 143                         for (unsigned ro=0; ro<size; ++ro) {
 144                                 // populate some elements
 145                                 ex elem = 0;
 146                                 if (rand()%(size/2) == 0)
 147                                         elem = sparse_tree(a, a, a, rand()%3, false, true, false);
 148                                 A.set(ro,co,elem);
 149                         }
 150                 }
 151                 ex det_gauss = A.determinant(determinant_algo::gauss);
 152                 ex det_laplace = A.determinant(determinant_algo::laplace);
 153                 ex det_divfree = A.determinant(determinant_algo::divfree);
 154                 ex det_bareiss = A.determinant(determinant_algo::bareiss);
 155                 if ((det_gauss-det_laplace).normal() != 0 ||
 156                         (det_bareiss-det_laplace).normal() != 0 ||
 157                         (det_divfree-det_laplace).normal() != 0) {
 158                         clog << "Determinant of " << size << "x" << size << " matrix "
 159                              << endl << A << endl
 160                              << "is inconsistent between different algorithms:" << endl
 161                              << "Gauss elimination:   " << det_gauss << endl
 162                              << "Minor elimination:   " << det_laplace << endl
 163                              << "Division-free elim.: " << det_divfree << endl
 164                              << "Fraction-free elim.: " << det_bareiss << endl;
 165                         ++result;
 166                 }
 167         }
 168
 169         return result;
 170 }
 171
 172 static unsigned symbolic_matrix_inverse()
 173 {
 174         unsigned result = 0;
 175         symbol a("a"), b("b"), c("c");
 176
 177         for (unsigned size=2; size<6; ++size) {
 178                 matrix A(size,size);
 179                 do {
 180                         for (unsigned co=0; co<size; ++co) {
 181                                 for (unsigned ro=0; ro<size; ++ro) {
 182                                         // populate some elements
 183                                         ex elem = 0;
 184                                         if (rand()%(size/2) == 0)
 185                                                 elem = sparse_tree(a, b, c, rand()%2, false, true, false);
 186                                         A.set(ro,co,elem);
 187                                 }
 188                         }
 189                 } while (A.determinant() == 0);
 190                 matrix B = A.inverse();
 191                 matrix C = A.mul(B);
 192                 bool ok = true;
 193                 for (unsigned ro=0; ro<size; ++ro)
 194                         for (unsigned co=0; co<size; ++co)
 195                                 if (C(ro,co).normal() != (ro==co?1:0))
 196                                         ok = false;
 197                 if (!ok) {
 198                         clog << "Inverse of " << size << "x" << size << " matrix "
 199                              << endl << A << endl
 200                              << "erroneously returned: "
 201                              << endl << B << endl;
 202                         ++result;
 203                 }
 204         }
 205
 206         return result;
 207 }
 208
 209 unsigned check_matrices()
 210 {
 211         unsigned result = 0;
 212
 213         cout << "checking symbolic matrix manipulations" << flush;
 214
 215         result += integdom_matrix_determinants();  cout << '.' << flush;
 216         result += rational_matrix_determinants();  cout << '.' << flush;
 217         result += funny_matrix_determinants();  cout << '.' << flush;
 218         result += compare_matrix_determinants();  cout << '.' << flush;
 219         result += symbolic_matrix_inverse();  cout << '.' << flush;
 220
 221         return result;
 222 }
 223
 224 int main(int argc, char** argv)
 225 {
 226         return check_matrices();
 227 }