self_avoiding_polygon/src/matrix.hpp

#ifndef MATRIX_HPP
#define MATRIX_HPP

#include <iostream>
#include <immintrin.h>

using namespace std;

class Matrix{
protected:
  size_t nr,nc,nc_avx,nc_full;
  double* data;
public:
  Matrix();
  void init(size_t n);
  void init(size_t nrow,size_t ncol);
  double get(size_t i,size_t j) const;
  double& get(size_t i,size_t j);
  void clear();
  void display() const;
  __m256d* get_avx_row(size_t i);
  const __m256d* get_avx_row(size_t i) const;

  void swap_lines(size_t i,size_t j);
  void mul_line(size_t i,double a);
  void add_mul_line(size_t i,size_t j,double a);
  double get_diag_square_sym(size_t i) const;
  double Gauss();
};

//*****************
//* AVX constants *
//*****************

//! Array of 4 zeros

static const __m256d zeros=_mm256_set1_pd(0);

//*************
//* Avx block *
//*************

//! Block structure mixing representing m256 as array of 4 doubles
union  AvxBlock{
  __m256d avx;
 double data[4];
};

//******************
//* Inline methods *
//******************

inline
Matrix::Matrix(){
  nr=0;
  nc=0;
  nc_avx=0;
  nc_full=0;
  data=nullptr;
}

inline void
Matrix::init(size_t n){
  return init(n,n);
}

inline double
Matrix::get(size_t i,size_t j) const{
  return data[i*nc_full+j];
}

inline double&
Matrix::get(size_t i,size_t j){
  return data[i*nc_full+j];
}

inline __m256d*
Matrix::get_avx_row(size_t r){
  return (__m256d*)(&data[r*nc_full]);
}

inline const __m256d*
Matrix::get_avx_row(size_t r) const{
  return (__m256d*)(&data[r*nc_full]);
}

#endif