// DENSEMat.H
#ifndef DENSEMAT_H
#define DENSEMAT_H
#include <string.h>
#include <cstdio>
#include <cstdlib>
#include "array.h"
#include "../dgtd-performance.hpp"

template <typename T>
class denseMat {
private:


public:
  int N, M;// N:row   M:column
  T *entry;
  explicit denseMat(int = 0, int = 0);
   denseMat(const denseMat& denMat);
  ~denseMat();
  // set functions
  void setDim(int, int);

  void addEntry(int rr, int cc, T value){ 
    entry[rr * M + cc] = entry[rr * M + cc] + value;
  }

  void setEntry(int rr, int cc, T value){
    entry[rr * M + cc] = value;
  }

  void reset(){ 
    for (int i = 0; i < N * M; i ++) 
      entry[i] = 0.0;
  }

  // get functions
  int getRowDim() const {return N;}
  int getColDim() const {return M;}
  T getEntry(int rr, int cc) const {return entry[rr * M + cc];}
  T* getEntryPtr() {return entry;}

  // mathematical operations BLAS 0 Operations
  denseMat& operator=(const denseMat &);
  denseMat operator+(const denseMat &) const;
  denseMat operator-(const denseMat &) const;
  denseMat operator*(const denseMat &) const;
  denseMat& operator*(const T &);
  // Vector operator*(const Vector& x) const;
  denseMat& operator*=(const T &);
  denseMat& operator+=(const denseMat &);
  denseMat& operator-=(const denseMat &);
  denseMat transpose();
  denseMat inverse();
  void SquareRoot(denseMat& SqRt);
  T FindSpectralRadius();
  void Clear();
  void scale(const T& val);
  void inverse(denseMat& invert);
  void transpose(denseMat* );
  void SelfTranspose();
  // print functions
  void print();
};
template class denseMat<fp_t>;

void solveAx_B(denseMat<fp_t>& A, ArrayFP<fp_t>& B); // solve Ax = B

void MXV(const denseMat<fp_t>, fp_t *, fp_t *);
void MXV(denseMat<fp_t>*, fp_t *, fp_t *);
void AXPY(denseMat<fp_t>* a, fp_t *x, fp_t *y);

// free functions
void Gaussj(denseMat<fp_t> *);
void Gaussj(denseMat<fp_t> *, bool&);

// BLAS 2 Operations
void aAxpby(fp_t alpha, fp_t beta, denseMat<fp_t>& A, fp_t* x, fp_t* y);
void aAxpby(fp_t alpha, fp_t beta, denseMat<fp_t>& A, ArrayFP<fp_t>& x, ArrayFP<fp_t>& y);
void aAxpy(fp_t alpha, denseMat<fp_t>& A, fp_t* x, fp_t* y);
void aAxpy(fp_t alpha, denseMat<fp_t>& A, ArrayFP<fp_t>& x, ArrayFP<fp_t>& y);
void Axpy(denseMat<fp_t>& A, fp_t* x, fp_t* y);
void Axpy(denseMat<fp_t>& A, ArrayFP<fp_t>& x, ArrayFP<fp_t>& y);
void aAx(fp_t alpha, denseMat<fp_t>& A, fp_t* x, fp_t* y);
void aAx(fp_t alpha, denseMat<fp_t>& A, ArrayFP<fp_t>& x, ArrayFP<fp_t>& y);
void Ax(denseMat<fp_t>& A, fp_t* x, fp_t* y);
void Ax(denseMat<fp_t>& A, ArrayFP<fp_t>& x, ArrayFP<fp_t>& y);

// BLAS 3 Operations
void aABpbC(fp_t alpha, fp_t beta, denseMat<fp_t>& A, denseMat<fp_t>& B, denseMat<fp_t>& C);
void aAB(fp_t alpha, denseMat<fp_t>& A, denseMat<fp_t>& B, denseMat<fp_t>& C);
void AB(denseMat<fp_t>& A, denseMat<fp_t>& B, denseMat<fp_t>& C);

#endif