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

#include "matrix.hpp"

void
Matrix::init(size_t nrow,size_t ncol){
  if(data!=nullptr){
    for(size_t i=0;i<nr*nc;++i) mpfr_clear(data[i]);
    delete[] data;
  }
  nr=nrow;
  nc=ncol;
  data=new mpfr_t[nr*nc];
  for(size_t i=0;i<nr*nc;++i) mpfr_init2(data[i],prec);
}

Matrix::~Matrix(){
  if(data!=nullptr){
    for(size_t i=0;i<nr*nc;++i) mpfr_clear(data[i]);
    delete[] data;
  }
}
//---------------
// Matrix::clear
//---------------

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

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

void
Matrix::mul_line(size_t i,mpfr_t x){
  for(size_t k=0;k<nc;++k){
    mpfr_mul(data[i*nc+k],data[i*nc+k],x,MPFR_RNDN);
  }
}

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

void
Matrix::Gauss(mpfr_t det){
  mpfr_set_ui(det,1,MPFR_RNDN);
  mpfr_t c;
  mpfr_init2(c,prec);
  size_t np=0; //np=0
  for(size_t j=0;j<nc;++j){
    for(size_t p=np;p<nr;++p){
      mpfr_set(c,get(p,j),MPFR_RNDN);
      if(mpfr_zero_p(c)==0){
        mpfr_mul(det,det,c,MPFR_RNDN);
        mpfr_ui_div(c,1,c,MPFR_RNDN);
        mul_line(p,c);
        for(size_t k=0;k<nr;++k){
          if(k!=p){
            mpfr_neg(c,get(k,j),MPFR_RNDN);
            add_mul_line(k,p,c);
          }
        }
        if(p!=np){
          swap_lines(np,p);
          mpfr_neg(det,det,MPFR_RNDN);
        }
        ++np;
        break;
      }
    }
  }
}