-/** Solve a set of equations for an m x n matrix.
- *
- * @param vars n x p matrix
- * @param rhs m x p matrix
- * @exception logic_error (incompatible matrices)
- * @exception runtime_error (singular matrix) */
-matrix matrix::solve(const matrix & vars,
- const matrix & rhs) const
-{
- if ((row != rhs.row) || (col != vars.row) || (rhs.col != vars.col))
- throw (std::logic_error("matrix::solve(): incompatible matrices"));
-
- throw (std::runtime_error("FIXME: need implementation."));
-}
-
-/** Old and obsolete interface: */
-matrix matrix::old_solve(const matrix & v) const
-{
- if ((v.row != col) || (col != v.row))
- throw (std::logic_error("matrix::solve(): incompatible matrices"));
-
- // build the augmented matrix of *this with v attached to the right
- matrix tmp(row,col+v.col);
- for (unsigned r=0; r<row; ++r) {
- for (unsigned c=0; c<col; ++c)
- tmp.m[r*tmp.col+c] = this->m[r*col+c];
- for (unsigned c=0; c<v.col; ++c)
- tmp.m[r*tmp.col+c+col] = v.m[r*v.col+c];
- }
- // cout << "augmented: " << tmp << endl;
- tmp.gauss_elimination();
- // cout << "degaussed: " << tmp << endl;
- // assemble the solution matrix
- exvector sol(v.row*v.col);
- for (unsigned c=0; c<v.col; ++c) {
- for (unsigned r=row; r>0; --r) {
- for (unsigned i=r; i<col; ++i)
- sol[(r-1)*v.col+c] -= tmp.m[(r-1)*tmp.col+i]*sol[i*v.col+c];
- sol[(r-1)*v.col+c] += tmp.m[(r-1)*tmp.col+col+c];
- sol[(r-1)*v.col+c] = (sol[(r-1)*v.col+c]/tmp.m[(r-1)*tmp.col+(r-1)]).normal();
- }
- }
- return matrix(v.row, v.col, sol);
-}
-