+#include <cstdlib> // for rand()
+#include <iostream>
+#include <sstream>
+using namespace std;
+
+extern const ex
+dense_univariate_poly(const symbol & x, unsigned degree);
+
+static unsigned check_matrix_solve(unsigned m, unsigned n, unsigned p,
+ unsigned degree)
+{
+ const symbol a("a");
+ matrix A(m,n);
+ matrix B(m,p);
+ // set the first min(m,n) rows of A and B
+ for (unsigned ro=0; (ro<m)&&(ro<n); ++ro) {
+ for (unsigned co=0; co<n; ++co)
+ A.set(ro,co,dense_univariate_poly(a,degree));
+ for (unsigned co=0; co<p; ++co)
+ B.set(ro,co,dense_univariate_poly(a,degree));
+ }
+ // repeat excessive rows of A and B to avoid excessive construction of
+ // overdetermined linear systems
+ for (unsigned ro=n; ro<m; ++ro) {
+ for (unsigned co=0; co<n; ++co)
+ A.set(ro,co,A(ro-1,co));
+ for (unsigned co=0; co<p; ++co)
+ B.set(ro,co,B(ro-1,co));
+ }
+ // create a vector of n*p symbols all named "xrc" where r and c are ints
+ vector<symbol> x;
+ matrix X(n,p);
+ for (unsigned i=0; i<n; ++i) {
+ for (unsigned j=0; j<p; ++j) {
+ ostringstream buf;
+ buf << "x" << i << j << ends;
+ x.push_back(symbol(buf.str()));
+ X.set(i,j,x[p*i+j]);
+ }
+ }
+ matrix sol(n,p);
+ // Solve the system A*X==B:
+ try {
+ sol = A.solve(X, B);
+ } catch (const exception & err) { // catch runtime_error
+ // Presumably, the coefficient matrix A was degenerate
+ string errwhat = err.what();
+ if (errwhat == "matrix::solve(): inconsistent linear system")
+ return 0;
+ else
+ clog << "caught exception: " << errwhat << endl;
+ throw;
+ }
+
+ // check the result with our original matrix:
+ bool errorflag = false;
+ for (unsigned ro=0; ro<m; ++ro) {
+ for (unsigned pco=0; pco<p; ++pco) {
+ ex e = 0;
+ for (unsigned co=0; co<n; ++co)
+ e += A(ro,co)*sol(co,pco);
+ if (!(e-B(ro,pco)).normal().is_zero())
+ errorflag = true;
+ }
+ }
+ if (errorflag) {
+ clog << "Our solve method claims that A*X==B, with matrices" << endl
+ << "A == " << A << endl
+ << "X == " << sol << endl
+ << "B == " << B << endl;
+ return 1;
+ }
+
+ return 0;
+}
+
+static unsigned check_inifcns_lsolve(unsigned n)