check/check_matrices.cpp

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