// DD_SPARMAT.H
#ifndef DD_SPARMAT_H
#define DD_SPARMAT_H
#include "densemat.h"
#include "hb_io.h"
#include "../dgtd-performance.hpp"

typedef enum {
  SYMMETRIC = 0,  // default
  ANTI_SYMMETRIC,
  NON_SYMMETRIC,
  HERMITIAN,
  SKEW_HERMITIAN
} MatrixType;

class sparMat {
  public:
    MatrixType type;
    int        dim;      // row dimension
    int        colDim;   // column dimension
    int        NZ;       // number of non zeros
    int        SEGSIZE;  // size of the segment
    int        NUMSEG;   // number of
    int*       rowA;     // rowA and colA are the Boeing format
    int*       colA;     // to store sparse matrix
    fp_t**   entry;    // at this moment, always entry[0], needs to be changed

  public:
    sparMat();
    ~sparMat();

    void init(int n, int nonZero);
    void sparInit(int, int);      // initialization of symmetric matrix
    void sparInit(int, int, int); // initialization of non-symmetric matrix
    void sparINIT(int, int, int); // initialization for non-symmetric matrix
    void sparInitNonSymm(int n, int nonZero);

    void memAlloc() { entry = new fp_t*[1]; entry[0] = new fp_t[NZ]; }
    void freeMemory();

  // modifiers
    void setMatType(MatrixType t)     { type = t; }
    void setDim(int n)                { dim = n; colDim = n; };
    void setRowDim(int n)             { dim = n; };
    void setColDim(int n)             { colDim = n; };
    void setNumSeg(int n)             { NUMSEG = n; };
    void setIA(int n, int ID);
    void setJA(int n, int ID);
    void setIAPtr(int *IAPTR)         { rowA = IAPTR; }
    void setJAPtr(int *JAPTR)         { colA = JAPTR; }
    void setEntry(int n, fp_t val) { entry[0][n] = val; }

  // accessors
    MatrixType isSymmetric()          { return type; }
    MatrixType isSymmetric()    const { return type; }
    MatrixType getMatType()           { return type; };
    MatrixType getMatType()     const { return type; }
    int*       getIAPtr()             { return rowA; }
    int*       getIAPtr()       const { return rowA; }
    int*       getJAPtr()             { return colA; }
    int*       getJAPtr()       const { return colA; }
    int        getRowMap(int i)       { return rowA[i]; }
    int        getRowMap(int i) const { return rowA[i]; }
    int        getColMap(int i)       { return colA[i]; };
    int        getColMap(int i) const { return colA[i]; };
    fp_t     getEntry(int n)        { return entry[0][n]; }
    fp_t     getEntry(int n)  const { return entry[0][n]; }
    fp_t*    getEntryPtr()          { return entry[0]; };
    fp_t*    getEntryPtr()    const { return entry[0]; };
    int        getDim()               { return dim; }
    int        getDim()         const { return dim; }
    int        getNZ()                { return NZ; }
    int        getNZ()          const { return NZ; }
    int        getRowDim()            { return dim; }
    int        getRowDim()      const { return dim; }
    int        getColDim()            { return colDim; }
    int        getColDim()      const { return colDim; }

  // utilities
    void addEntry(int row, int col, fp_t val);
    void blockAddEntry(int rDim, int cDim, int rowMap[], int colMap[], fp_t **entry);
    void blockAddEntryNonSymm(int rDim, int cDim, int rowMap[], int colMap[], denseMat<fp_t> &D);
    void blockAddEntry(int rDim, int cDim, int rowMap[], int colMap[], denseMat<fp_t> &D);
    void blockAddEntry(int nn, int map[], denseMat<fp_t>& D);
    void blockAddEntrySymm(int rDim, int cDim, int rowMap[], int colMap[], const denseMat<fp_t> &D);
    void blockAddEntryAnti(int rDim, int cDim, int rowMap[], int colMap[], const denseMat<fp_t> &D);
    void diagonalScale(fp_t *d1, fp_t *d2);
    sparMat* symmToFull();

  // I/O
    void print();
    void print(char *fname);
    void printIA(char *fname);
    void printJA(char *fname);
    void preRead(char *fname);
    void read(char *fname);
    void hb_file_write(char const* fname);
    void bmtsave(const char*, int base = 0, int float_prec  =64,int int_prec = 32) const;
    int  loadBmtSparse(const char *fname, char *order, int *m, int *n, int *nnz, fp_t **val,
                  int **ind, int **ptr);
    void st_to_hb(int ,int, int ,int [],int [],fp_t []);

    fp_t memUsage();
    //
    void NonSymmMxV(fp_t* ,fp_t* );
};

// /////////////////////////////////////////
// // free functions
// Vector MxV(sparMat *const A, Vector &x);
// Vector MxV(sparMat const &A, Vector &x);
// Vector MxV(sparMat *const A, Vector &x, int *const xMap, int *const yMap);
// 
// Vector MxV(sparMat *const A, Vector &x, bool transpose);
// void symmMxV(sparMat *const A, Vector &x, Vector &y);
// void nonSymmMxV(sparMat *const A, Vector &x,  Vector &y, bool transpose);
   
// void symm_MxV(sparMat const &A, Vector &x, Vector &y);
// void nonsymm_MxV(sparMat const &A, Vector &x, Vector &y);
// void antisymm_MxV(sparMat const &A, Vector &x, Vector &y);
// void hermitian_MxV(sparMat const &A, Vector &x, Vector &y);
// void skew_herm_MxV(sparMat const &A, Vector &x, Vector &y);
// 
// Vector operator*(const sparMat &A, Vector &x);

#endif