self_avoiding_polygon/other/trash/float128/src/matrix.cpp

#include "matrix.hpp"

void
Matrix::init(size_t nrow,size_t ncol){
  if(data!=nullptr) delete[] data;
  nr=nrow;
  nc=ncol;
  data=new Reel[nr*nc];
}

//---------------
// Matrix::clear
//---------------

void
Matrix::clear(){
  for(size_t i=0;i<nr*nc;++i){
    data[i]=0;
  }
}

void
Matrix::display() const{
  /*for(size_t i=0;i<nr;++i){
    for(size_t j=0;j<nc;++j){
      cout<<get(i,j)<<'\t';
    }
    cout<<endl;
  }*/
}

void
Matrix::swap_lines(size_t i,size_t j){
  for(size_t k=0;k<nc;++k){
    swap(data[i*nc+k],data[j*nc+k]);
  }
}

void
Matrix::mul_line(size_t i,Reel a){
  for(size_t k=0;k<nc;++k){
    data[i*nc+k]*=a;
  }
}

void
Matrix::add_mul_line(size_t i,size_t j,Reel a){
  for(size_t k=0;k<nc;++k){
    data[i*nc+k]+=a*data[j*nc+k];
  }
}


Reel
Matrix::Gauss(){
  Reel det=1;
  size_t np=0; //np=0
  for(size_t j=0;j<nc;++j){
    for(size_t p=np;p<nr;++p){
      Reel c=get(p,j);
      if(c!=0){
        det*=c;
        mul_line(p,1.0/c);
        for(size_t k=0;k<nr;++k){
          if(k!=p){
            add_mul_line(k,p,-get(k,j));
          }
        }
        if(p!=np){
          swap_lines(np,p);
          det*=-1;
        }
        ++np;
        break;
      }
    }
  }
  return det;
}

Reel
Matrix::get_diag_square_sym(size_t i) const{
  Reel res=0;
  for(size_t k=0;k<nc;++k){
    res+=get(i,k);
  }
  return res;
}