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