SMatrix.h

#ifndef S_MATRIX_H
#define S_MATRIX_H
 
#include <tuple>
#include <vector>
 
#include "Blas.h"
#include "exceptions.h"
 
template <typename T>
class SerialMatrix {
  int numRows_;
  int numCols_;
  size_t numElements_;
 
  T* mat = nullptr;  // pointer to the internal array structure.
 
  std::tuple<int, int> local2Global(const int& k);
 
 public:
  int global2Local(const int& row, const int& col);
 
  static const char transN = 'N';
  static const char transT = 'T';
  static const char transC = 'C';
 
  SerialMatrix(const int& numRows, const int& numCols,
                 const int& numBlasRows = 0, const int& numBlasCols = 0,
                 const int& numBlocksRows = 0, const int& numBlocksCols = 0);
 
  SerialMatrix();
 
  ~SerialMatrix();
 
  SerialMatrix(const SerialMatrix<T>& that);
 
  SerialMatrix<T>& operator=(const SerialMatrix<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;
  size_t size() const;
  T& operator()(const int &row, const int &col);
  const T& operator()(const int &row, const int &col) const;
 
  void enforcePositiveSemiDefinite() {}; 
 
  void outputToHDF5(const std::string &outFileName,
                                   const std::string &dataSetName) {}
 
  SerialMatrix<T> prod(const SerialMatrix<T>& that, const char& trans1 = transN,
                 const char& trans2 = transN);
 
  SerialMatrix<T> operator+=(const SerialMatrix<T>& m1) {
    if(numRows_ != m1.rows() || numCols_ != m1.cols()) {
      Error("Cannot add matrices of different sizes.");
    }
    for (size_t s = 0; s < size(); s++) mat[s] += m1.mat[s];
    return *this;
  }
 
  SerialMatrix<T> operator-=(const SerialMatrix<T>& m1) {
    if(numRows_ != m1.rows() || numCols_ != m1.cols()) {
      Error("Cannot subtract matrices of different sizes.");
    }
    for (size_t s = 0; s < size(); s++) mat[s] -= m1.mat[s];
    return *this;
  }
 
  SerialMatrix<T> operator*=(const T& that) {
    for (size_t s = 0; s < size(); s++) mat[s] *= that;
    return *this;
  }
 
  SerialMatrix<T> operator/=(const T& that) {
    for (size_t s = 0; s < size(); s++) mat[s] /= that;
    return *this;
  }
 
  void eye();
 
  std::tuple<std::vector<double>, SerialMatrix<T>> diagonalize();
 
  std::tuple<std::vector<double>, SerialMatrix<T>> diagonalize(int numEigenvalues,
                                                bool checkNegativeEigenvalues = true);
 
  double squaredNorm();
 
  double norm();
 
  T dot(const SerialMatrix<T>& that);
 
  SerialMatrix<T> operator-() const;
 
  void symmetrize();
 
};
 
// A default constructor to build a dense matrix of zeros to be filled
template <typename T>
SerialMatrix<T>::SerialMatrix(const int& numRows, const int& numCols,
                  const int& numBlasRows, const int& numBlasCols,
                  const int& numBlocksRows, const int& numBlocksCols) {
  // these are used only the in the pmatrix case
  (void) numBlocksRows;
  (void) numBlocksCols;
  (void) numBlasRows;
  (void) numBlasCols;
 
  numRows_ = numRows;
  numCols_ = numCols;
  numElements_ = numRows_ * numCols_;
  mat = new T[numRows_ * numCols_];
  for (size_t i = 0; i < numElements_; i++) mat[i] = 0;  // fill with zeroes
  assert(mat != nullptr);  // Memory could not be allocated, end program
}
 
// default constructor
template <typename T>
SerialMatrix<T>::SerialMatrix() {
  mat = nullptr;
  numRows_ = 0;
  numCols_ = 0;
  numElements_ = 0;
}
 
// copy constructor
template <typename T>
SerialMatrix<T>::SerialMatrix(const SerialMatrix<T>& that) {
  numRows_ = that.numRows_;
  numCols_ = that.numCols_;
  numElements_ = that.numElements_;
  if (mat != nullptr) {
    delete[] mat;
    mat = nullptr;
  }
  mat = new T[numElements_];
  assert(mat != nullptr);
  for (size_t i = 0; i < numElements_; i++) {
    mat[i] = that.mat[i];
  }
}
 
template <typename T>
SerialMatrix<T>& SerialMatrix<T>::operator=(const SerialMatrix<T>& that) {
  if (this != &that) {
    numRows_ = that.numRows_;
    numCols_ = that.numCols_;
    numElements_ = that.numElements_;
    // matrix allocation
    if (mat != nullptr) {
      delete[] mat;
    }
    mat = new T[numElements_];
    assert(mat != nullptr);
    for (size_t i = 0; i < numElements_; i++) {
      mat[i] = that.mat[i];
    }
  }
  return *this;
}
 
// destructor
template <typename T>
SerialMatrix<T>::~SerialMatrix() {
  if (mat != nullptr) delete[] mat;
}
 
/* ------------- Very basic operations -------------- */
template <typename T>
int SerialMatrix<T>::rows() const {
  return numRows_;
}
template <typename T>
int SerialMatrix<T>::localRows() const {
  return numRows_;
}
template <typename T>
int SerialMatrix<T>::cols() const {
  return numCols_;
}
template <typename T>
int SerialMatrix<T>::localCols() const {
  return numCols_;
}
template <typename T>
size_t SerialMatrix<T>::size() const {
  return numElements_;
}
 
// Get/set element
template <typename T>
T& SerialMatrix<T>::operator()(const int &row, const int &col) {
  if(row >= numRows_ || col >= numCols_ || row < 0 || col < 0) {
    DeveloperError("Tried to reference a SMatrix state out of bounds: " 
        + std::to_string(row) + " " + std::to_string(col));
  }
  return mat[global2Local(row, col)];
}
 
template <typename T>
const T& SerialMatrix<T>::operator()(const int &row, const int &col) const {
  if(row >= numRows_ || col >= numCols_ || row < 0 || col < 0) {
    DeveloperError("Tried to reference a SMatrix state out of bounds: " 
        + std::to_string(row) + " " + std::to_string(col));
  }
  return mat[global2Local(row, col)];
}
 
template <typename T>
bool SerialMatrix<T>::indicesAreLocal(const int& row, const int& col) {
  if(row >= numRows_ || col >= numCols_ || row < 0 || col < 0) {
    DeveloperError("indicesAreLocal tried to reference a SMatrix state out of bounds: " 
          + std::to_string(row) + " " + std::to_string(col));
  }
  return true;
}
 
template <typename T>
std::tuple<int, int> SerialMatrix<T>::local2Global(const int& k) {
  // we convert this combined local index k into row / col indices
  // k = j * nRows + i
  if(numRows_ == 0) Error("attempted to div by zero in l2g");
  int j = k / numRows_;
  int i = k - j * numRows_;
  return std::make_tuple(i, j);
}
 
// Indexing to set up the matrix in col major format
template <typename T>
int SerialMatrix<T>::global2Local(const int& row, const int& col) {
  return numRows_ * col + row;
}
 
template <typename T>
std::vector<std::tuple<int, int>> SerialMatrix<T>::getAllLocalStates() {
  std::vector<std::tuple<int, int>> x;
  for (size_t k = 0; k < numElements_; k++) {
    std::tuple<int, int> t = local2Global(k);  // bloch indices
    x.push_back(t);
  }
  return x;
}
 
// General unary negation
template <typename T>
SerialMatrix<T> SerialMatrix<T>::operator-() const {
  SerialMatrix<T> ret(numRows_, numCols_);
  for (int row = 0; row < numRows_; row++) {
    for (int col = 0; col < numCols_; col++) ret(row, col) = -(*this)(row, col);
  }
  return ret;
}
 
// Sets the matrix to the identity matrix
template <typename T>
void SerialMatrix<T>::eye() {
  if(numRows_ != numCols_) {
    Error("Cannot build an identity matrix with non-square matrix");
  }
  for (int row = 0; row < numRows_; row++) (*this)(row, row) = (T)1.0;
}
 
// forward declaration of explicit specialization
// is needed because the generic implementation does not
// work for complex types
template<> double SerialMatrix<std::complex<double>>::norm();
 
template <typename T>
double SerialMatrix<T>::norm() {
  T sumSq = 0;
  for (int row = 0; row < numRows_; row++) {
    for (int col = 0; col < numCols_; col++) {
      sumSq += ((*this)(row, col) * (*this)(row, col));
    }
  }
  return sqrt(sumSq);
}
 
template <typename T>
double SerialMatrix<T>::squaredNorm() {
  double x = norm();
  return x * x;
}
 
template <typename T>
T SerialMatrix<T>::dot(const SerialMatrix<T>& that) {
  T scalar = (T)0.;
  for (size_t i = 0; i < numElements_; i++) {
    scalar += (*(mat + i)) * (*(that.mat + i));
  }
  return scalar;
}
 
#endif  // S_MATRIX_H