check/check_lsolve.cpp

   1 /** @file check_lsolve.cpp
   2  *
   3  *  These test routines do some simple checks on solving linear systems of
   4  *  symbolic equations. */
   5
   6 /*
   7  *  GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany
   8  *
   9  *  This program is free software; you can redistribute it and/or modify
  10  *  it under the terms of the GNU General Public License as published by
  11  *  the Free Software Foundation; either version 2 of the License, or
  12  *  (at your option) any later version.
  13  *
  14  *  This program is distributed in the hope that it will be useful,
  15  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  17  *  GNU General Public License for more details.
  18  *
  19  *  You should have received a copy of the GNU General Public License
  20  *  along with this program; if not, write to the Free Software
  21  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  22  */
  23
  24 #include "checks.h"
  25
  26 static unsigned lsolve1(int size)
  27 {
  28     // A dense size x size matrix in dense univariate random polynomials
  29     // of order 4.
  30     unsigned result = 0;
  31     symbol a("a");
  32     ex sol;
  33
  34     // Create two dense linear matrices A and B where all entries are random
  35     // univariate polynomials
  36     matrix A(size,size), B(size,2), X(size,2);
  37     for (int ro=0; ro<size; ++ro) {
  38         for (int co=0; co<size; ++co)
  39             A.set(ro,co,dense_univariate_poly(a, 5));
  40         for (int co=0; co<2; ++co)
  41             B.set(ro,co,dense_univariate_poly(a, 5));
  42     }
  43     if (A.determinant().is_zero())
  44         clog << "lsolve1: singular system!" << endl;
  45
  46     // Solve the system A*X==B:
  47     X = A.old_solve(B);
  48
  49     // check the result:
  50     bool errorflag = false;
  51     matrix Aux(size,2);
  52     Aux = A.mul(X).sub(B);
  53     for (int ro=0; ro<size && !errorflag; ++ro)
  54         for (int co=0; co<2; ++co)
  55             if (!(Aux(ro,co)).normal().is_zero())
  56                 errorflag = true;
  57     if (errorflag) {
  58         clog << "Our solve method claims that A*X==B, with matrices" << endl
  59              << "A == " << A << endl
  60              << "X == " << X << endl
  61              << "B == " << B << endl;
  62         ++result;
  63     }
  64     return result;
  65 }
  66
  67 static unsigned lsolve2(int size)
  68 {
  69     // A dense size x size matrix in dense bivariate random polynomials
  70     // of order 2.
  71     unsigned result = 0;
  72     symbol a("a"), b("b");
  73     ex sol;
  74
  75     // Create two dense linear matrices A and B where all entries are dense random
  76     // bivariate polynomials:
  77     matrix A(size,size), B(size,2), X(size,2);
  78     for (int ro=0; ro<size; ++ro) {
  79         for (int co=0; co<size; ++co)
  80             A.set(ro,co,dense_bivariate_poly(a,b,2));
  81         for (int co=0; co<2; ++co)
  82             B.set(ro,co,dense_bivariate_poly(a,b,2));
  83     }
  84     if (A.determinant().is_zero())
  85         clog << "lsolve2: singular system!" << endl;
  86
  87     // Solve the system A*X==B:
  88     X = A.old_solve(B);
  89
  90     // check the result:
  91     bool errorflag = false;
  92     matrix Aux(size,2);
  93     Aux = A.mul(X).sub(B);
  94     for (int ro=0; ro<size && !errorflag; ++ro)
  95         for (int co=0; co<2; ++co)
  96             if (!(Aux(ro,co)).normal().is_zero())
  97                 errorflag = true;
  98     if (errorflag) {
  99         clog << "Our solve method claims that A*X==B, with matrices" << endl
 100              << "A == " << A << endl
 101              << "X == " << X << endl
 102              << "B == " << B << endl;
 103         ++result;
 104     }
 105     return result;
 106 }
 107
 108 unsigned check_lsolve(void)
 109 {
 110     unsigned result = 0;
 111
 112     cout << "checking linear solve" << flush;
 113     clog << "---------linear solve:" << endl;
 114
 115     //result += lsolve1(2);  cout << '.' << flush;
 116     //result += lsolve1(3);  cout << '.' << flush;
 117     //result += lsolve2(2);  cout << '.' << flush;
 118     //result += lsolve2(3);  cout << '.' << flush;
 119
 120     if (!result) {
 121         cout << " passed " << endl;
 122         clog << "(no output)" << endl;
 123     } else {
 124         cout << " failed " << endl;
 125     }
 126
 127     return result;
 128 }