Matrix.h

#ifndef MATRIX_H
#define MATRIX_H
 
#include "PMatrix.h"
#include "SMatrix.h"
 
template <typename T>
class Matrix {
 
  int global2Local(const int& row, const int& col);
  std::tuple<int, int> local2Global(const int& k);
 
 bool isDistributed;
 
 ParallelMatrix<T>* pmat = nullptr;
 
 SerialMatrix<T>* mat = nullptr;
 
 public:
  Matrix(const int& numRows, const int& numCols, const int& numBlocksRows = 0,
         const int& numBlocksCols = 0, bool isDistributed_ = false);
 
  Matrix();
 
  ~Matrix();
 
  Matrix(const Matrix<T>& that);
 
  Matrix<T>& operator=(const Matrix<T>& that);
 
  std::vector<std::tuple<int, int>> getAllLocalStates();
 
  bool indicesAreLocal(const int& row, const int& col);
 
  int rows() const;
  int localRows() const;
  int cols() const;
  int localCols() const;
  int size() const;
  double getMemory() const;
 
  T& operator()(const int &row, const int &col);
 
  const T& operator()(const int &row, const int &col) const;
 
  Matrix<T> prod(const Matrix<T>& that, const char& trans1,
                 const char& trans2);
 
  Matrix<T> operator+=(const Matrix<T>& m1) {
    if(isDistributed) (*pmat) += (*m1.pmat);
    else (*mat) += (*m1.mat);
    return *this;
  }
 
  Matrix<T> operator-=(const Matrix<T>& m1) {
    if(isDistributed) (*pmat) -= (*m1.pmat);
    else (*mat) -= (*m1.mat);
    return *this;
  }
 
  Matrix<T> operator*=(const T& that) {
    if(isDistributed) (*pmat) *= that;
    else (*mat) *= that;
    return *this;
  }
 
  Matrix<T> operator/=(const T& that) {
    if(isDistributed) (*pmat) /= that;
    else (*mat) /= that;
    return *this;
  }
 
  void eye();
 
  std::tuple<std::vector<double>, Matrix<T>> diagonalize();
 
  std::tuple<std::vector<double>, Matrix<T>> diagonalize(int numEigenvalues,
                                                bool checkNegativeEigenvalues = true);
 
  double squaredNorm();
 
  double norm();
 
  T dot(const Matrix<T>& that);
 
  Matrix<T> operator-() const;
 
  void outputToHDF5();
 
  void symmetrize();
 
};
 
/* ------------------ constructor implementations -------------- */
 
// A default constructor to build a dense matrix of zeros to be filled
template <typename T>
Matrix<T>::Matrix(const int& numRows, const int& numCols,
                  const int& numBlocksRows, const int& numBlocksCols, bool isDistributed_) {
 
  isDistributed = isDistributed_; // default to false if no value supplied
 
  if(isDistributed){
    pmat = new ParallelMatrix<T>(numRows,numCols,numBlocksRows,numBlocksCols);
  }
  else {
    mat = new SerialMatrix<T>(numRows,numCols);
  }
}
 
// default constructor
template <typename T>
Matrix<T>::Matrix() {
  isDistributed = false;
  if (pmat!=nullptr) delete pmat;
  if (mat!=nullptr) delete mat;
  mat = new SerialMatrix<T>();
}
 
// copy constructor
template <typename T>
Matrix<T>::Matrix(const Matrix<T>& that) {
  isDistributed = that.isDistributed;
 
  // call SMatrix or PMatrix copy constructor
  if(isDistributed) {
    pmat = that.pmat;
  }
  else {
    mat = that.mat;
  }
}
 
template <typename T>
Matrix<T>& Matrix<T>::operator=(const Matrix<T>& that) {
  if (this != &that) {
    isDistributed = that.isDistributed;
    // call SMatrix or PMatrix copy constructor
    if(isDistributed) {
      if (pmat!=nullptr) delete pmat;
      pmat = new ParallelMatrix<T>(*that.pmat);
    }
    else {
      if (mat!=nullptr) delete mat;
      mat = new SerialMatrix<T>(*that.mat);
    }
  }
  return *this;
}
 
// destructor
template <typename T>
Matrix<T>::~Matrix() {
  if (pmat!=nullptr) {
    delete pmat;
  }
  if (mat!=nullptr) {
    delete mat;
  }
}
 
/* ------------- Very basic operations -------------- */
template <typename T>
int Matrix<T>::rows() const {
  if(isDistributed) return pmat->rows();
  else{ return mat->rows(); }
}
template <typename T>
int Matrix<T>::localRows() const {
  if(isDistributed) return pmat->localRows();
  else{ return mat->rows(); }
}
template <typename T> int Matrix<T>::cols() const {
  if(isDistributed) return pmat->cols();
  else{ return mat->cols(); }
}
template <typename T>
int Matrix<T>::localCols() const {
  if(isDistributed) return pmat->localCols();
  else{ return mat->cols(); }
}
template <typename T>
int Matrix<T>::size() const {
  if(isDistributed) return pmat->size();
  else{ return mat->size(); }
}
 
/* ------------- get-set operations -------------- */
template <typename T>
double Matrix<T>::getMemory() const{
  // this is done in parts to avoid overflow;
  // size in GB is size of type*rows()*cols()/(1024)^3
  double temp = sizeof(T)*rows()/1024.;
  temp = temp*cols()/1024.;
  return temp/(1024.);
}
 
template <typename T>
T& Matrix<T>::operator()(const int &row, const int &col) {
  if(isDistributed) return (*pmat)(row,col);
  else { return (*mat)(row,col); }
}
 
template <typename T>
const T& Matrix<T>::operator()(const int &row, const int &col) const {
  if(isDistributed) return (*pmat)(row,col);
  else{ return (*mat)(row,col); }
}
 
template <typename T>
bool Matrix<T>::indicesAreLocal(const int& row, const int& col) {
  if(isDistributed) return pmat->indicesAreLocal(row,col);
  else{ return true; }
}
 
template <typename T>
std::tuple<int, int> Matrix<T>::local2Global(const int& k) {
  if(isDistributed) return pmat->local2Global(k);
  else{ return mat->local2Global(k); }
}
 
// Indexing to set up the matrix in col major format
template <typename T>
int Matrix<T>::global2Local(const int& row, const int& col) {
  if(isDistributed) return pmat->global2Local(row,col);
  else{ return mat->global2Local(row,col); }
}
 
template <typename T>
std::vector<std::tuple<int, int>> Matrix<T>::getAllLocalStates() {
  if(isDistributed) return pmat->getAllLocalStates();
  else{ return mat->getAllLocalStates(); }
}
 
/* ------------- basic linear algebra ops -------------- */
// General unary negation
template <typename T>
Matrix<T> Matrix<T>::operator-() const {
  Matrix<T> c(*this); // copy this matrix
    if (isDistributed) *(c.pmat) = -*(c.pmat);
    else               *(c.mat)  = -*(c.mat);
  return c;
}
 
// Sets the matrix to the identity matrix
template <typename T>
void Matrix<T>::eye() {
  if(isDistributed) pmat->eye();
  else{ mat->eye(); }
}
 
template <typename T>
double Matrix<T>::norm() {
  if(isDistributed) return pmat->norm();
  else{ return mat->norm(); }
}
 
template <typename T>
double Matrix<T>::squaredNorm() {
  if(isDistributed) return pmat->squaredNorm();
  else{ return mat->squaredNorm(); }
}
 
template <typename T>
T Matrix<T>::dot(const Matrix<T>& that) {
  if(isDistributed) return pmat->dot(that.pmat);
  else{ return mat->dot(that.mat); }
}
 
template <typename T>
void Matrix<T>::outputToHDF5() {
  if(isDistributed) pmat->outputToHDF5();
  else{ Error("Write to HDF5 not implemented for SMatrix."); }
}
 
#endif  // MATRIX_H