- 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) {
- char buf[4];
- ostrstream(buf,sizeof(buf)) << i << j << ends;
- x.push_back(symbol(string("x")+buf));
- 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;
+ 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;