- // A dense size x size matrix in dense univariate random polynomials
- // of order 4.
- unsigned result = 0;
- symbol a("a");
- ex sol;
-
- // Create two dense linear matrices A and B where all entries are random
- // univariate polynomials
- matrix A(size,size), B(size,2), X(size,2);
- for (int ro=0; ro<size; ++ro) {
- for (int co=0; co<size; ++co)
- A.set(ro,co,dense_univariate_poly(a, 5));
- for (int co=0; co<2; ++co)
- B.set(ro,co,dense_univariate_poly(a, 5));
- }
- if (A.determinant().is_zero())
- clog << "lsolve1: singular system!" << endl;
-
- // Solve the system A*X==B:
- X = A.old_solve(B);
-
- // check the result:
- bool errorflag = false;
- matrix Aux(size,2);
- Aux = A.mul(X).sub(B);
- for (int ro=0; ro<size && !errorflag; ++ro)
- for (int co=0; co<2; ++co)
- if (!(Aux(ro,co)).normal().is_zero())
- errorflag = true;
- if (errorflag) {
- clog << "Our solve method claims that A*X==B, with matrices" << endl
- << "A == " << A << endl
- << "X == " << X << endl
- << "B == " << B << endl;
- ++result;
- }
- return result;
+ 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;