ginac/matrix.cpp

   1 /** @file matrix.cpp
   2  *
   3  *  Implementation of symbolic matrices
   4  *
   5  *  GiNaC Copyright (C) 1999 Johannes Gutenberg University Mainz, Germany
   6  *
   7  *  This program is free software; you can redistribute it and/or modify
   8  *  it under the terms of the GNU General Public License as published by
   9  *  the Free Software Foundation; either version 2 of the License, or
  10  *  (at your option) any later version.
  11  *
  12  *  This program is distributed in the hope that it will be useful,
  13  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15  *  GNU General Public License for more details.
  16  *
  17  *  You should have received a copy of the GNU General Public License
  18  *  along with this program; if not, write to the Free Software
  19  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  20  */
  21
  22 #include <algorithm>
  23 #include <stdexcept>
  24
  25 #include "matrix.h"
  26
  27 //////////
  28 // default constructor, destructor, copy constructor, assignment operator
  29 // and helpers:
  30 //////////
  31
  32 // public
  33
  34 /** Default ctor.  Initializes to 1 x 1-dimensional zero-matrix. */
  35 matrix::matrix()
  36     : basic(TINFO_matrix), row(1), col(1)
  37 {
  38     debugmsg("matrix default constructor",LOGLEVEL_CONSTRUCT);
  39     m.push_back(exZERO());
  40 }
  41
  42 matrix::~matrix()
  43 {
  44     debugmsg("matrix destructor",LOGLEVEL_DESTRUCT);
  45 }
  46
  47 matrix::matrix(matrix const & other)
  48 {
  49     debugmsg("matrix copy constructor",LOGLEVEL_CONSTRUCT);
  50     copy(other);
  51 }
  52
  53 matrix const & matrix::operator=(matrix const & other)
  54 {
  55     debugmsg("matrix operator=",LOGLEVEL_ASSIGNMENT);
  56     if (this != &other) {
  57         destroy(1);
  58         copy(other);
  59     }
  60     return *this;
  61 }
  62
  63 // protected
  64
  65 void matrix::copy(matrix const & other)
  66 {
  67     basic::copy(other);
  68     row=other.row;
  69     col=other.col;
  70     m=other.m;  // use STL's vector copying
  71 }
  72
  73 void matrix::destroy(bool call_parent)
  74 {
  75     if (call_parent) basic::destroy(call_parent);
  76 }
  77
  78 //////////
  79 // other constructors
  80 //////////
  81
  82 // public
  83
  84 /** Very common ctor.  Initializes to r x c-dimensional zero-matrix.
  85  *
  86  *  @param r number of rows
  87  *  @param c number of cols */
  88 matrix::matrix(int r, int c)
  89     : basic(TINFO_matrix), row(r), col(c)
  90 {
  91     debugmsg("matrix constructor from int,int",LOGLEVEL_CONSTRUCT);
  92     m.resize(r*c, exZERO());
  93 }
  94
  95 // protected
  96
  97 /** Ctor from representation, for internal use only. */
  98 matrix::matrix(int r, int c, vector<ex> const & m2)
  99     : basic(TINFO_matrix), row(r), col(c), m(m2)
 100 {
 101     debugmsg("matrix constructor from int,int,vector<ex>",LOGLEVEL_CONSTRUCT);
 102 }
 103
 104 //////////
 105 // functions overriding virtual functions from bases classes
 106 //////////
 107
 108 // public
 109
 110 basic * matrix::duplicate() const
 111 {
 112     debugmsg("matrix duplicate",LOGLEVEL_DUPLICATE);
 113     return new matrix(*this);
 114 }
 115
 116 /** nops is defined to be rows x columns. */
 117 int matrix::nops() const
 118 {
 119     return row*col;
 120 }
 121
 122 /** returns matrix entry at position (i/col, i%col). */
 123 ex & matrix::let_op(int const i)
 124 {
 125     return m[i];
 126 }
 127
 128 /** expands the elements of a matrix entry by entry. */
 129 ex matrix::expand(unsigned options) const
 130 {
 131     vector<ex> tmp(row*col);
 132     for (int i=0; i<row*col; ++i) {
 133         tmp[i]=m[i].expand(options);
 134     }
 135     return matrix(row, col, tmp);
 136 }
 137
 138 /** Search ocurrences.  A matrix 'has' an expression if it is the expression
 139  *  itself or one of the elements 'has' it. */
 140 bool matrix::has(ex const & other) const
 141 {
 142     ASSERT(other.bp!=0);
 143
 144     // tautology: it is the expression itself
 145     if (is_equal(*other.bp)) return true;
 146
 147     // search all the elements
 148     for (vector<ex>::const_iterator r=m.begin(); r!=m.end(); ++r) {
 149         if ((*r).has(other)) return true;
 150     }
 151     return false;
 152 }
 153
 154 /** evaluate matrix entry by entry. */
 155 ex matrix::eval(int level) const
 156 {
 157     debugmsg("matrix eval",LOGLEVEL_MEMBER_FUNCTION);
 158
 159     // check if we have to do anything at all
 160     if ((level==1)&&(flags & status_flags::evaluated)) {
 161         return *this;
 162     }
 163
 164     // emergency break
 165     if (level == -max_recursion_level) {
 166         throw (std::runtime_error("matrix::eval(): recursion limit exceeded"));
 167     }
 168
 169     // eval() entry by entry
 170     vector<ex> m2(row*col);
 171     --level;
 172     for (int r=0; r<row; ++r) {
 173         for (int c=0; c<col; ++c) {
 174             m2[r*col+c] = m[r*col+c].eval(level);
 175         }
 176     }
 177
 178     return (new matrix(row, col, m2))->setflag(status_flags::dynallocated |
 179                                                status_flags::evaluated );
 180 }
 181
 182 /** evaluate matrix numerically entry by entry. */
 183 ex matrix::evalf(int level) const
 184 {
 185     debugmsg("matrix evalf",LOGLEVEL_MEMBER_FUNCTION);
 186
 187     // check if we have to do anything at all
 188     if (level==1) {
 189         return *this;
 190     }
 191
 192     // emergency break
 193     if (level == -max_recursion_level) {
 194         throw (std::runtime_error("matrix::evalf(): recursion limit exceeded"));
 195     }
 196
 197     // evalf() entry by entry
 198     vector<ex> m2(row*col);
 199     --level;
 200     for (int r=0; r<row; ++r) {
 201         for (int c=0; c<col; ++c) {
 202             m2[r*col+c] = m[r*col+c].evalf(level);
 203         }
 204     }
 205     return matrix(row, col, m2);
 206 }
 207
 208 // protected
 209
 210 int matrix::compare_same_type(basic const & other) const
 211 {
 212     ASSERT(is_exactly_of_type(other, matrix));
 213     matrix const & o=static_cast<matrix &>(const_cast<basic &>(other));
 214
 215     // compare number of rows
 216     if (row != o.rows()) {
 217         return row < o.rows() ? -1 : 1;
 218     }
 219
 220     // compare number of columns
 221     if (col != o.cols()) {
 222         return col < o.cols() ? -1 : 1;
 223     }
 224
 225     // equal number of rows and columns, compare individual elements
 226     int cmpval;
 227     for (int r=0; r<row; ++r) {
 228         for (int c=0; c<col; ++c) {
 229             cmpval=((*this)(r,c)).compare(o(r,c));
 230             if (cmpval!=0) return cmpval;
 231         }
 232     }
 233     // all elements are equal => matrices are equal;
 234     return 0;
 235 }
 236
 237 //////////
 238 // non-virtual functions in this class
 239 //////////
 240
 241 // public
 242
 243 /** Sum of matrices.
 244  *
 245  *  @exception logic_error (incompatible matrices) */
 246 matrix matrix::add(matrix const & other) const
 247 {
 248     if (col != other.col || row != other.row) {
 249         throw (std::logic_error("matrix::add(): incompatible matrices"));
 250     }
 251
 252     vector<ex> sum(this->m);
 253     vector<ex>::iterator i;
 254     vector<ex>::const_iterator ci;
 255     for (i=sum.begin(), ci=other.m.begin();
 256          i!=sum.end();
 257          ++i, ++ci) {
 258         (*i) += (*ci);
 259     }
 260     return matrix(row,col,sum);
 261 }
 262
 263 /** Difference of matrices.
 264  *
 265  *  @exception logic_error (incompatible matrices) */
 266 matrix matrix::sub(matrix const & other) const
 267 {
 268     if (col != other.col || row != other.row) {
 269         throw (std::logic_error("matrix::sub(): incompatible matrices"));
 270     }
 271
 272     vector<ex> dif(this->m);
 273     vector<ex>::iterator i;
 274     vector<ex>::const_iterator ci;
 275     for (i=dif.begin(), ci=other.m.begin();
 276          i!=dif.end();
 277          ++i, ++ci) {
 278         (*i) -= (*ci);
 279     }
 280     return matrix(row,col,dif);
 281 }
 282
 283 /** Product of matrices.
 284  *
 285  *  @exception logic_error (incompatible matrices) */
 286 matrix matrix::mul(matrix const & other) const
 287 {
 288     if (col != other.row) {
 289         throw (std::logic_error("matrix::mul(): incompatible matrices"));
 290     }
 291
 292     vector<ex> prod(row*other.col);
 293     for (int i=0; i<row; ++i) {
 294         for (int j=0; j<other.col; ++j) {
 295             for (int l=0; l<col; ++l) {
 296                 prod[i*other.col+j] += m[i*col+l] * other.m[l*other.col+j];
 297             }
 298         }
 299     }
 300     return matrix(row, other.col, prod);
 301 }
 302
 303 /** operator() to access elements.
 304  *
 305  *  @param ro row of element
 306  *  @param co column of element
 307  *  @exception range_error (index out of range) */
 308 ex const & matrix::operator() (int ro, int co) const
 309 {
 310     if (ro<0 || ro>=row || co<0 || co>=col) {
 311         throw (std::range_error("matrix::operator(): index out of range"));
 312     }
 313
 314     return m[ro*col+co];
 315 }
 316
 317 /** Set individual elements manually.
 318  *
 319  *  @exception range_error (index out of range) */
 320 matrix & matrix::set(int ro, int co, ex value)
 321 {
 322     if (ro<0 || ro>=row || co<0 || co>=col) {
 323         throw (std::range_error("matrix::set(): index out of range"));
 324     }
 325
 326     ensure_if_modifiable();
 327     m[ro*col+co]=value;
 328     return *this;
 329 }
 330
 331 /** Transposed of an m x n matrix, producing a new n x m matrix object that
 332  *  represents the transposed. */
 333 matrix matrix::transpose(void) const
 334 {
 335     vector<ex> trans(col*row);
 336
 337     for (int r=0; r<col; ++r) {
 338         for (int c=0; c<row; ++c) {
 339             trans[r*row+c] = m[c*col+r];
 340         }
 341     }
 342     return matrix(col,row,trans);
 343 }
 344
 345 /* Determiant of purely numeric matrix, using pivoting. This routine is only
 346  * called internally by matrix::determinant(). */
 347 ex determinant_numeric(const matrix & M)
 348 {
 349     ASSERT(M.rows()==M.cols());  // cannot happen, just in case...
 350     matrix tmp(M);
 351     ex det=exONE();
 352     ex piv;
 353
 354     for (int r1=0; r1<M.rows(); ++r1) {
 355         int indx = tmp.pivot(r1);
 356         if (indx == -1) {
 357             return exZERO();
 358         }
 359         if (indx != 0) {
 360             det *= exMINUSONE();
 361         }
 362         det = det * tmp.m[r1*M.cols()+r1];
 363         for (int r2=r1+1; r2<M.rows(); ++r2) {
 364             piv = tmp.m[r2*M.cols()+r1] / tmp.m[r1*M.cols()+r1];
 365             for (int c=r1+1; c<M.cols(); c++) {
 366                 tmp.m[r2*M.cols()+c] -= piv * tmp.m[r1*M.cols()+c];
 367             }
 368         }
 369     }
 370     return det;
 371 }
 372
 373 // Compute the sign of a permutation of a vector of things, used internally
 374 // by determinant_symbolic_perm() where it is instantiated for int.
 375 template <class T>
 376 int permutation_sign(vector<T> s)
 377 {
 378     if (s.size() < 2)
 379         return 0;
 380     int sigma=1;
 381     for (typename vector<T>::iterator i=s.begin(); i!=s.end()-1; ++i) {
 382         for (typename vector<T>::iterator j=i+1; j!=s.end(); ++j) {
 383             if (*i == *j)
 384                 return 0;
 385             if (*i > *j) {
 386                 iter_swap(i,j);
 387                 sigma = -sigma;
 388             }
 389         }
 390     }
 391     return sigma;
 392 }
 393
 394 /** Determinant built by application of the full permutation group. This
 395  *  routine is only called internally by matrix::determinant(). */
 396 ex determinant_symbolic_perm(const matrix & M)
 397 {
 398     ASSERT(M.rows()==M.cols());  // cannot happen, just in case...
 399
 400     if (M.rows()==1) {  // speed things up
 401         return M(0,0);
 402     }
 403
 404     ex det;
 405     ex term;
 406     vector<int> sigma(M.cols());
 407     for (int i=0; i<M.cols(); ++i) sigma[i]=i;
 408
 409     do {
 410         term = M(sigma[0],0);
 411         for (int i=1; i<M.cols(); ++i) term *= M(sigma[i],i);
 412         det += permutation_sign(sigma)*term;
 413     } while (next_permutation(sigma.begin(), sigma.end()));
 414
 415     return det;
 416 }
 417
 418 /** Recursive determiant for small matrices having at least one symbolic entry.
 419  *  This algorithm is also known as Laplace-expansion. This routine is only
 420  *  called internally by matrix::determinant(). */
 421 ex determinant_symbolic_minor(const matrix & M)
 422 {
 423     ASSERT(M.rows()==M.cols());  // cannot happen, just in case...
 424
 425     if (M.rows()==1) {  // end of recursion
 426         return M(0,0);
 427     }
 428     if (M.rows()==2) {  // speed things up
 429         return (M(0,0)*M(1,1)-
 430                 M(1,0)*M(0,1));
 431     }
 432     if (M.rows()==3) {  // speed things up even a little more
 433         return ((M(2,1)*M(0,2)-M(2,2)*M(0,1))*M(1,0)+
 434                 (M(1,2)*M(0,1)-M(1,1)*M(0,2))*M(2,0)+
 435                 (M(2,2)*M(1,1)-M(2,1)*M(1,2))*M(0,0));
 436     }
 437
 438     ex det;
 439     matrix minorM(M.rows()-1,M.cols()-1);
 440     for (int r1=0; r1<M.rows(); ++r1) {
 441         // assemble the minor matrix
 442         for (int r=0; r<minorM.rows(); ++r) {
 443             for (int c=0; c<minorM.cols(); ++c) {
 444                 if (r<r1) {
 445                     minorM.set(r,c,M(r,c+1));
 446                 } else {
 447                     minorM.set(r,c,M(r+1,c+1));
 448                 }
 449             }
 450         }
 451         // recurse down
 452         if (r1%2) {
 453             det -= M(r1,0) * determinant_symbolic_minor(minorM);
 454         } else {
 455             det += M(r1,0) * determinant_symbolic_minor(minorM);
 456         }
 457     }
 458     return det;
 459 }
 460
 461 /*  Leverrier algorithm for large matrices having at least one symbolic entry.
 462  *  This routine is only called internally by matrix::determinant(). The
 463  *  algorithm is deemed bad for symbolic matrices since it returns expressions
 464  *  that are very hard to canonicalize. */
 465 /*ex determinant_symbolic_leverrier(const matrix & M)
 466  *{
 467  *    ASSERT(M.rows()==M.cols());  // cannot happen, just in case...
 468  *
 469  *    matrix B(M);
 470  *    matrix I(M.row, M.col);
 471  *    ex c=B.trace();
 472  *    for (int i=1; i<M.row; ++i) {
 473  *        for (int j=0; j<M.row; ++j)
 474  *            I.m[j*M.col+j] = c;
 475  *        B = M.mul(B.sub(I));
 476  *        c = B.trace()/ex(i+1);
 477  *    }
 478  *    if (M.row%2) {
 479  *        return c;
 480  *    } else {
 481  *        return -c;
 482  *    }
 483  *}*/
 484
 485 /** Determinant of square matrix.  This routine doesn't actually calculate the
 486  *  determinant, it only implements some heuristics about which algorithm to
 487  *  call.  When the parameter for normalization is explicitly turned off this
 488  *  method does not normalize its result at the end, which might imply that
 489  *  the symbolic 2x2 matrix [[a/(a-b),1],[b/(a-b),1]] is not immediatly
 490  *  recognized to be unity.  (This is Mathematica's default behaviour, it
 491  *  should be used with care.)
 492  *
 493  *  @param     normalized may be set to false if no normalization of the
 494  *             result is desired (i.e. to force Mathematica behavior, Maple
 495  *             does normalize the result).
 496  *  @return    the determinant as a new expression
 497  *  @exception logic_error (matrix not square) */
 498 ex matrix::determinant(bool normalized) const
 499 {
 500     if (row != col) {
 501         throw (std::logic_error("matrix::determinant(): matrix not square"));
 502     }
 503
 504     // check, if there are non-numeric entries in the matrix:
 505     for (vector<ex>::const_iterator r=m.begin(); r!=m.end(); ++r) {
 506         if (!(*r).info(info_flags::numeric)) {
 507             if (normalized) {
 508                 return determinant_symbolic_minor(*this).normal();
 509             } else {
 510                 return determinant_symbolic_perm(*this);
 511             }
 512         }
 513     }
 514     // if it turns out that all elements are numeric
 515     return determinant_numeric(*this);
 516 }
 517
 518 /** Trace of a matrix.
 519  *
 520  *  @return    the sum of diagonal elements
 521  *  @exception logic_error (matrix not square) */
 522 ex matrix::trace(void) const
 523 {
 524     if (row != col) {
 525         throw (std::logic_error("matrix::trace(): matrix not square"));
 526     }
 527
 528     ex tr;
 529     for (int r=0; r<col; ++r) {
 530         tr += m[r*col+r];
 531     }
 532     return tr;
 533 }
 534
 535 /** Characteristic Polynomial.  The characteristic polynomial of a matrix M is
 536  *  defined as the determiant of (M - lambda * 1) where 1 stands for the unit
 537  *  matrix of the same dimension as M.  This method returns the characteristic
 538  *  polynomial as a new expression.
 539  *
 540  *  @return    characteristic polynomial as new expression
 541  *  @exception logic_error (matrix not square)
 542  *  @see       matrix::determinant() */
 543 ex matrix::charpoly(ex const & lambda) const
 544 {
 545     if (row != col) {
 546         throw (std::logic_error("matrix::charpoly(): matrix not square"));
 547     }
 548
 549     matrix M(*this);
 550     for (int r=0; r<col; ++r) {
 551         M.m[r*col+r] -= lambda;
 552     }
 553     return (M.determinant());
 554 }
 555
 556 /** Inverse of this matrix.
 557  *
 558  *  @return    the inverted matrix
 559  *  @exception logic_error (matrix not square)
 560  *  @exception runtime_error (singular matrix) */
 561 matrix matrix::inverse(void) const
 562 {
 563     if (row != col) {
 564         throw (std::logic_error("matrix::inverse(): matrix not square"));
 565     }
 566
 567     matrix tmp(row,col);
 568     // set tmp to the unit matrix
 569     for (int i=0; i<col; ++i) {
 570         tmp.m[i*col+i] = exONE();
 571     }
 572     // create a copy of this matrix
 573     matrix cpy(*this);
 574     for (int r1=0; r1<row; ++r1) {
 575         int indx = cpy.pivot(r1);
 576         if (indx == -1) {
 577             throw (std::runtime_error("matrix::inverse(): singular matrix"));
 578         }
 579         if (indx != 0) {  // swap rows r and indx of matrix tmp
 580             for (int i=0; i<col; ++i) {
 581                 tmp.m[r1*col+i].swap(tmp.m[indx*col+i]);
 582             }
 583         }
 584         ex a1 = cpy.m[r1*col+r1];
 585         for (int c=0; c<col; ++c) {
 586             cpy.m[r1*col+c] /= a1;
 587             tmp.m[r1*col+c] /= a1;
 588         }
 589         for (int r2=0; r2<row; ++r2) {
 590             if (r2 != r1) {
 591                 ex a2 = cpy.m[r2*col+r1];
 592                 for (int c=0; c<col; ++c) {
 593                     cpy.m[r2*col+c] -= a2 * cpy.m[r1*col+c];
 594                     tmp.m[r2*col+c] -= a2 * tmp.m[r1*col+c];
 595                 }
 596             }
 597         }
 598     }
 599     return tmp;
 600 }
 601
 602 void matrix::ffe_swap(int r1, int c1, int r2 ,int c2)
 603 {
 604     ensure_if_modifiable();
 605
 606     ex tmp=ffe_get(r1,c1);
 607     ffe_set(r1,c1,ffe_get(r2,c2));
 608     ffe_set(r2,c2,tmp);
 609 }
 610
 611 void matrix::ffe_set(int r, int c, ex e)
 612 {
 613     set(r-1,c-1,e);
 614 }
 615
 616 ex matrix::ffe_get(int r, int c) const
 617 {
 618     return operator()(r-1,c-1);
 619 }
 620
 621 /** Solve a set of equations for an m x n matrix by fraction-free Gaussian
 622  *  elimination. Based on algorithm 9.1 from 'Algorithms for Computer Algebra'
 623  *  by Keith O. Geddes et al.
 624  *
 625  *  @param vars n x p matrix
 626  *  @param rhs m x p matrix
 627  *  @exception logic_error (incompatible matrices)
 628  *  @exception runtime_error (singular matrix) */
 629 matrix matrix::fraction_free_elim(matrix const & vars,
 630                                   matrix const & rhs) const
 631 {
 632     if ((row != rhs.row) || (col != vars.row) || (rhs.col != vars.col)) {
 633         throw (std::logic_error("matrix::solve(): incompatible matrices"));
 634     }
 635
 636     matrix a(*this); // make a copy of the matrix
 637     matrix b(rhs);     // make a copy of the rhs vector
 638
 639     // given an m x n matrix a, reduce it to upper echelon form
 640     int m=a.row;
 641     int n=a.col;
 642     int sign=1;
 643     ex divisor=1;
 644     int r=1;
 645
 646     // eliminate below row r, with pivot in column k
 647     for (int k=1; (k<=n)&&(r<=m); ++k) {
 648         // find a nonzero pivot
 649         int p;
 650         for (p=r; (p<=m)&&(a.ffe_get(p,k).is_equal(exZERO())); ++p) {}
 651         // pivot is in row p
 652         if (p<=m) {
 653             if (p!=r) {
 654                 // switch rows p and r
 655                 for (int j=k; j<=n; ++j) {
 656                     a.ffe_swap(p,j,r,j);
 657                 }
 658                 b.ffe_swap(p,1,r,1);
 659                 // keep track of sign changes due to row exchange
 660                 sign=-sign;
 661             }
 662             for (int i=r+1; i<=m; ++i) {
 663                 for (int j=k+1; j<=n; ++j) {
 664                     a.ffe_set(i,j,(a.ffe_get(r,k)*a.ffe_get(i,j)
 665                                   -a.ffe_get(r,j)*a.ffe_get(i,k))/divisor);
 666                     a.ffe_set(i,j,a.ffe_get(i,j).normal() /*.normal() */ );
 667                 }
 668                 b.ffe_set(i,1,(a.ffe_get(r,k)*b.ffe_get(i,1)
 669                               -b.ffe_get(r,1)*a.ffe_get(i,k))/divisor);
 670                 b.ffe_set(i,1,b.ffe_get(i,1).normal() /*.normal() */ );
 671                 a.ffe_set(i,k,0);
 672             }
 673             divisor=a.ffe_get(r,k);
 674             r++;
 675         }
 676     }
 677     // optionally compute the determinant for square or augmented matrices
 678     // if (r==m+1) { det=sign*divisor; } else { det=0; }
 679
 680     /*
 681     for (int r=1; r<=m; ++r) {
 682         for (int c=1; c<=n; ++c) {
 683             cout << a.ffe_get(r,c) << "\t";
 684         }
 685         cout << " | " <<  b.ffe_get(r,1) << endl;
 686     }
 687     */
 688
 689 #ifdef DOASSERT
 690     // test if we really have an upper echelon matrix
 691     int zero_in_last_row=-1;
 692     for (int r=1; r<=m; ++r) {
 693         int zero_in_this_row=0;
 694         for (int c=1; c<=n; ++c) {
 695             if (a.ffe_get(r,c).is_equal(exZERO())) {
 696                zero_in_this_row++;
 697             } else {
 698                 break;
 699             }
 700         }
 701         ASSERT((zero_in_this_row>zero_in_last_row)||(zero_in_this_row=n));
 702         zero_in_last_row=zero_in_this_row;
 703     }
 704 #endif // def DOASSERT
 705
 706     // assemble solution
 707     matrix sol(n,1);
 708     int last_assigned_sol=n+1;
 709     for (int r=m; r>0; --r) {
 710         int first_non_zero=1;
 711         while ((first_non_zero<=n)&&(a.ffe_get(r,first_non_zero).is_zero())) {
 712             first_non_zero++;
 713         }
 714         if (first_non_zero>n) {
 715             // row consists only of zeroes, corresponding rhs must be 0 as well
 716             if (!b.ffe_get(r,1).is_zero()) {
 717                 throw (std::runtime_error("matrix::fraction_free_elim(): singular matrix"));
 718             }
 719         } else {
 720             // assign solutions for vars between first_non_zero+1 and
 721             // last_assigned_sol-1: free parameters
 722             for (int c=first_non_zero+1; c<=last_assigned_sol-1; ++c) {
 723                 sol.ffe_set(c,1,vars.ffe_get(c,1));
 724             }
 725             ex e=b.ffe_get(r,1);
 726             for (int c=first_non_zero+1; c<=n; ++c) {
 727                 e=e-a.ffe_get(r,c)*sol.ffe_get(c,1);
 728             }
 729             sol.ffe_set(first_non_zero,1,
 730                         (e/a.ffe_get(r,first_non_zero)).normal());
 731             last_assigned_sol=first_non_zero;
 732         }
 733     }
 734     // assign solutions for vars between 1 and
 735     // last_assigned_sol-1: free parameters
 736     for (int c=1; c<=last_assigned_sol-1; ++c) {
 737         sol.ffe_set(c,1,vars.ffe_get(c,1));
 738     }
 739
 740     /*
 741     for (int c=1; c<=n; ++c) {
 742         cout << vars.ffe_get(c,1) << "->" << sol.ffe_get(c,1) << endl;
 743     }
 744     */
 745
 746 #ifdef DOASSERT
 747     // test solution with echelon matrix
 748     for (int r=1; r<=m; ++r) {
 749         ex e=0;
 750         for (int c=1; c<=n; ++c) {
 751             e=e+a.ffe_get(r,c)*sol.ffe_get(c,1);
 752         }
 753         if (!(e-b.ffe_get(r,1)).normal().is_zero()) {
 754             cout << "e=" << e;
 755             cout << "b.ffe_get(" << r<<",1)=" << b.ffe_get(r,1) << endl;
 756             cout << "diff=" << (e-b.ffe_get(r,1)).normal() << endl;
 757         }
 758         ASSERT((e-b.ffe_get(r,1)).normal().is_zero());
 759     }
 760
 761     // test solution with original matrix
 762     for (int r=1; r<=m; ++r) {
 763         ex e=0;
 764         for (int c=1; c<=n; ++c) {
 765             e=e+ffe_get(r,c)*sol.ffe_get(c,1);
 766         }
 767         try {
 768         if (!(e-rhs.ffe_get(r,1)).normal().is_zero()) {
 769             cout << "e=" << e << endl;
 770             e.printtree(cout);
 771             ex en=e.normal();
 772             cout << "e.normal()=" << en << endl;
 773             en.printtree(cout);
 774             cout << "rhs.ffe_get(" << r<<",1)=" << rhs.ffe_get(r,1) << endl;
 775             cout << "diff=" << (e-rhs.ffe_get(r,1)).normal() << endl;
 776         }
 777         } catch (...) {
 778             ex xxx=e-rhs.ffe_get(r,1);
 779             cerr << "xxx=" << xxx << endl << endl;
 780         }
 781         ASSERT((e-rhs.ffe_get(r,1)).normal().is_zero());
 782     }
 783 #endif // def DOASSERT
 784
 785     return sol;
 786 }
 787
 788 /** Solve simultaneous set of equations. */
 789 matrix matrix::solve(matrix const & v) const
 790 {
 791     if (!(row == col && col == v.row)) {
 792         throw (std::logic_error("matrix::solve(): incompatible matrices"));
 793     }
 794
 795     // build the extended matrix of *this with v attached to the right
 796     matrix tmp(row,col+v.col);
 797     for (int r=0; r<row; ++r) {
 798         for (int c=0; c<col; ++c) {
 799             tmp.m[r*tmp.col+c] = m[r*col+c];
 800         }
 801         for (int c=0; c<v.col; ++c) {
 802             tmp.m[r*tmp.col+c+col] = v.m[r*v.col+c];
 803         }
 804     }
 805     for (int r1=0; r1<row; ++r1) {
 806         int indx = tmp.pivot(r1);
 807         if (indx == -1) {
 808             throw (std::runtime_error("matrix::solve(): singular matrix"));
 809         }
 810         for (int c=r1; c<tmp.col; ++c) {
 811             tmp.m[r1*tmp.col+c] /= tmp.m[r1*tmp.col+r1];
 812         }
 813         for (int r2=r1+1; r2<row; ++r2) {
 814             for (int c=r1; c<tmp.col; ++c) {
 815                 tmp.m[r2*tmp.col+c]
 816                     -= tmp.m[r2*tmp.col+r1] * tmp.m[r1*tmp.col+c];
 817             }
 818         }
 819     }
 820
 821     // assemble the solution matrix
 822     vector<ex> sol(v.row*v.col);
 823     for (int c=0; c<v.col; ++c) {
 824         for (int r=col-1; r>=0; --r) {
 825             sol[r*v.col+c] = tmp[r*tmp.col+c];
 826             for (int i=r+1; i<col; ++i) {
 827                 sol[r*v.col+c]
 828                     -= tmp[r*tmp.col+i] * sol[i*v.col+c];
 829             }
 830         }
 831     }
 832     return matrix(v.row, v.col, sol);
 833 }
 834
 835 // protected
 836
 837 /** Partial pivoting method.
 838  *  Usual pivoting returns the index to the element with the largest absolute
 839  *  value and swaps the current row with the one where the element was found.
 840  *  Here it does the same with the first non-zero element. (This works fine,
 841  *  but may be far from optimal for numerics.) */
 842 int matrix::pivot(int ro)
 843 {
 844     int k=ro;
 845
 846     for (int r=ro; r<row; ++r) {
 847         if (!m[r*col+ro].is_zero()) {
 848             k = r;
 849             break;
 850         }
 851     }
 852     if (m[k*col+ro].is_zero()) {
 853         return -1;
 854     }
 855     if (k!=ro) {  // swap rows
 856         for (int c=0; c<col; ++c) {
 857             m[k*col+c].swap(m[ro*col+c]);
 858         }
 859         return k;
 860     }
 861     return 0;
 862 }
 863
 864 //////////
 865 // global constants
 866 //////////
 867
 868 const matrix some_matrix;
 869 type_info const & typeid_matrix=typeid(some_matrix);