4d1f23730aeda0acb20b388ec522a692835acb75
[ginac.git] / 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 }