matrix.H

Go to the documentation of this file.

#ifndef _MATRIX_H_
#define _MATRIX_H_

#include "tables.H"
using namespace std;


template <class T> class Matrix : public Table<T> {
public:
  Matrix<T>() : Table<T>() {this->left='('; this->right=')';}
  Matrix<T>(long i, long j) : Table<T>(i,j) {this->left='('; this->right=')';}
  Matrix<T>(const Matrix<T>& that) : Table<T>(that) {}

  void make_identity(long n)
  {
    this->resize(n, n);
    long i, j;
    for(i= 0; i< this->nb_lines; i++)
      for(j= 0; j< this->nb_cols; j++)
        if( i == j )
          this->pdata[i][i]= 1;
        else
          this->pdata[i][j]= 0;
  }


  void get_data_from_matrix_list(const list< Matrix<T> >& L){
    long li, col, i, j, total;
    typename list< Matrix<T> >::const_iterator it;
    
    //find out number of lines and cols
    li= col= 0;
    for(it= L.begin(); it!= L.end(); it++){
      li+= it->nlines();
      col= col > it->ncols() ? col : it->ncols();
    }
    //get the memory for it
    delete[] this->data;
    this->data= new T[col*li];
    this->nb_lines= li;
    this->nb_cols= col;
    this->tuple_size= col*li;
    delete[] this->pdata;
    this->write_ptrs();

    //now copy
    for(total= 0, it=L.begin(); it!= L.end(); total+= it->nlines(), it++)
      for(i=0; i< it->nlines(); i++)
        for(j=0; j< it->ncols(); j++)
          this->pdata[total+i][j]= (*it)(i,j);

  }

  

  void swap_lines(long i, long j){
    long k;
    T dummy;
    
    if( i == j )
      return;

    for(k=0; k< this->nb_cols; k++){
      dummy= this->pdata[i][k];
      this->pdata[i][k]=  this->pdata[j][k];
      this->pdata[j][k]= dummy;
    }
  }
  void add_line(long i, long j, T coeff){
    long k;

    for(k=0; k< this->nb_cols; k++)
      this->pdata[i][k]+= this->pdata[j][k]*coeff;
  }

  void multiply_line(long i, T coeff)
  {
    for(long j=0; j< this->nb_cols; j++)
      this->pdata[i][j]*= coeff;
  }
  
    


  bool is_invertible() const {
    Matrix<T> M;
    long col, line;
        

    if( this->nb_cols != this->nb_lines )
      return false;

    M= *this;
    for(col=0; col < this->nb_cols; col++){

      // look for a nonzero entry
      for(line= col; line< this->nb_lines; line++)
        if( M(line, col) != 0 )
          break;
      

      if(line == this->nb_lines) //not found ?
        return false;
      // else swap the lines
      M.swap_lines(line, col);
      // and triangulate
      for(line= col + 1; line < this->nb_lines; line++)
        M.add_line(line, col, - M(line, col) / M(col, col) );
    }
    return true;

  }

  Matrix<T> inverse() const 
  {
    if(this->nb_cols != this->nb_lines)
      return *this;
    //else...
    Matrix<T> M, N;
    M= *this;
    N.make_identity(this->nb_lines);
    long line, col;
    T c;
    
        
    for(col=0; col < this->nb_cols; col++){
      //looks for a nonzero entry in M
      for(line= col; line< this->nb_lines; line++)
        if( M(line, col) != 0 )
          break;
      if(line == this->nb_lines) //not found ?
        return *this;
      c= 1 / M(line,col);
      M.multiply_line(line, c);
      N.multiply_line(line, c);
      M.swap_lines(line, col);
      N.swap_lines(line, col);
      for(line=0; line<this->nb_lines; line++)
        if(line != col){
          c= -M(line, col);
          M.add_line(line, col, c);
          N.add_line(line, col, c);
        }
    
    }
    
    return N;
  }
  



  Matrix<T>& operator+=(const Matrix<T>& that){
    long i;
    
    if( (this->nb_lines != that.nb_lines) || (this->nb_cols != that.nb_cols) )
      return *this;
    for(i=0; i < this->tuple_size; i++){
      this->data[i] += that.data[i];
    }
    return *this;
  }

   Matrix<T>& operator-=(const Matrix<T>& that){
    long i,j;
    
    if( (this->nb_lines != that.nb_lines) || (this->nb_cols != that.nb_cols) )
      return *this;
    for(i=0; i < this->tuple_size; i++){
      this->data[i] -= that.data[i];
    }
    return *this;
  }



  Matrix<T>& operator*=(const Matrix<T>& that){
    Matrix<T> ans;
    
    matrixprod(*this, that, ans);
    this->steal( ans );
    return *this;
  }
    //OLD VERSION
//     Tuple<T> temp;
//     Matrix<T> *duplicate;
//     T element;
//     long current_line, current_col,j;

//     if( this->nb_cols != that.nb_lines )
//       return *this;
    
//     if( this == &that ){   // if we are squaring this...
//       duplicate=new Matrix<T>(*this);
//       (*this) *= (*duplicate);
//       delete duplicate;
//       return *this;
//     }

//     this->resize(this->nb_lines, that.nb_cols);
//     temp.resize(that.nb_cols);

//     for(current_line=0; current_line < this->nb_lines; current_line++){
//       for(current_col=0; current_col < that.nb_cols; current_col++){
//      element= this->pdata[current_line][0];
//      element *= that.pdata[0][current_col];
//      temp[current_col]= element;
//      for(j=1; j < that.nb_lines; j++){
//        element= this->pdata[current_line][j];
//        element *= that.pdata[j][current_col];
//        temp[current_col] += element;
//      }
//       }
//       readline(current_line, temp);
//     }

//     return *this;
//  }

  Matrix<T>& operator*=(const T& c) {
    long i;

    for(i=0; i< this->nb_lines * this->nb_cols; i++)
      this->data[i] *= c;
     
    return *this;
  }

  Matrix<T> operator*(const T& c) const {
    Matrix<T> ans;

    ans= *this;
    ans*= c;
    return ans;
  }




friend void matrixprod(const Matrix<T>& A, const Matrix<T>& B, Matrix<T>& ans) {
    T element;
    long current_line, current_col,j;

    if( A.nb_cols != B.nb_lines)
      return;
    
    ans.resize( A.nb_lines , B.nb_cols );

    for(current_line=0; current_line < A.nb_lines; current_line++){
      for(current_col=0; current_col< B.nb_cols; current_col++){
        element = A(current_line,0);
        element  *= B(0,current_col);   
        ans(current_line,current_col)= element;
        for(j=1; j < A.nb_cols; j++){
          element = A(current_line,j);
          element  *= B(j,current_col); 
          ans(current_line,current_col) += element;
        }
      }
    }

    return;
  }


  Matrix<T> operator*(const Matrix<T>& that) const {
    Matrix<T> ans;

    matrixprod(*this, that, ans);
    return ans;
  }

  Matrix<T> operator+(const Matrix<T>& that) const {
    Matrix<T> ans;

    ans= *this;
    ans+= that;
    return ans;
  }


};

#endif

---