Program Listing for File qmatrix.hpp

↰ Return to documentation for file (utils/qmatrix.hpp)

#pragma once
// #include "utils/timer.hpp"
#include <cmath>
#include <complex>
#include <initializer_list>
#include <iostream>
#include <numeric>
// #include <lapacke.h> // comment this if mkl is used
// #include <cblas.h>
// uncomment if you want to use the mkl class
#define MKL_Complex16 std::complex<double>
#include <mkl.h>
#include <mkl_lapacke.h>
// typedef double _Complex DCOMPLEX
// #include <omp.h>
#include <ostream>
#include <stdexcept>
#include <tuple>
#include <type_traits>
#include <vector>
template <typename NType1>
std::ostream &operator<<(std::ostream &out, const std::vector<NType1> &val) {
  for (auto const &aa : val) {
    out << aa << " ";
  }
  out << "\n";
  return out;
}
template <typename NType1>
std::ostream &operator<<(std::ostream                           &out,
                         const std::vector<std::vector<NType1>> &val) {
  out << "\n";
  for (auto x : val) {
    for (auto y : x) {
      out << y << " ";
    }
    out << "\n";
  }
  return out;
}
template <typename T> void LOGGER(const std::string &name, T x) {
  std::cout << "## " << name << " : " << x << std::endl;
}
// #define logger(name) LOGGER(#name, (name))
// This is a matrix class quamtum(q) operators
// This is also consistent with Lapack_row_MAJOR
using cm_vec = std::vector<std::complex<double>>;
template <class T = double> // default is double
class qmatrix {
  std::vector<T> mat{0};
  size_t         row{0}, column{0}, dim{0};

public:
  qmatrix(const std::initializer_list<T> inputVec, size_t _row = 0,
          size_t _column = 0)
      : mat(inputVec) {
    // check if this is a quare matrix
    // std::cout << "constructed with a " << inputVec.size() << "-element
    // list\n";
    if (_row == 0 && _column == 0) {
      this->row    = std::sqrt(inputVec.size());
      this->column = std::sqrt(inputVec.size());
      this->dim    = inputVec.size();
      if (this->dim != this->row * this->column) {
        throw std::runtime_error(
            "initializer_list: vector is not a square matrix! ");
      }
    } else {
      this->row    = _row;
      this->column = _column;
      this->dim    = _column * _row;
    }
    this->mat.shrink_to_fit();
  }
  qmatrix(const std::vector<T> &inputVec, size_t _row = 0, // NOLINT
          size_t _column = 0) {
    // check if this is a quare matrix
    if (_row == 0 && _column == 0) {
      this->row    = std::sqrt(inputVec.size());
      this->column = std::sqrt(inputVec.size());
      this->dim    = inputVec.size();
      if (this->dim != this->row * this->column) {
        throw std::runtime_error("vector is not a square matrix! ");
      }
    } else {
      this->row    = _row;
      this->column = _column;
      this->dim    = _column * _row;
    }
    this->mat = inputVec;
    this->mat.shrink_to_fit();
  }
  qmatrix(size_t _row, size_t _column, T populate) {
    this->row    = _row;
    this->column = _column;
    this->dim    = _column * _row;
    this->mat    = std::vector<T>(_row * _column, populate);
    this->mat.shrink_to_fit();
    // TODO(sp): restrict memory usage
  }
  void resize(size_t _row = 0, size_t _column = 0, T populate = 0) {
    this->row    = _row;
    this->column = _column;
    this->dim    = _column * _row;
    this->mat    = std::vector<T>(_row * _column, populate);
    this->mat.shrink_to_fit();
    // TODO(sp): restrict memory usage
  }
  void clear() {
    this->row    = 0;
    this->column = 0;
    this->dim    = 0;
    this->mat.clear();
    this->mat.shrink_to_fit();
  } // This is for square matrix
  qmatrix(size_t N, T populate) { qmatrix(N, N, populate); }
  qmatrix() { qmatrix(0, 0, 0); }
  [[nodiscard]] T     &operator()(size_t i) { return this->mat[i]; }
  [[nodiscard]] T      operator()(size_t i) const { return this->mat[i]; }
  [[nodiscard]] T     &at(size_t i) { return this->mat[i]; }
  [[nodiscard]] T      at(size_t i) const { return this->mat[i]; }
  [[nodiscard]] size_t size() const { return dim; }
  [[nodiscard]] size_t getrow() const { return row; }
  [[nodiscard]] size_t getcolumn() const { return column; }
  [[nodiscard]] T     &operator()(size_t i, size_t j) {
    return this->mat[i * column + j];
  }
  [[nodiscard]] T operator()(size_t i, size_t j) const {
    return this->mat[i * column + j];
  }
  // add a at operator
  [[nodiscard]] T &at(size_t i, size_t j) { return this->mat[i * column + j]; }
  [[nodiscard]] T  at(size_t i, size_t j) const {
    return this->mat[i * column + j];
  }
  // TODO(sp): arithmetic operator
  // TODO(sp): use stl
  [[nodiscard]] T sum() const {
    // Sum of all elements
    return std::accumulate(this->mat.begin(), this->mat.end(), 0);
  }
  [[nodiscard]] T absSum() const {
    // Sum of all absolute values of the elements
    double sum2 = 0;
    for (auto aa : this->mat) {
      sum2 = sum2 + std::fabs(aa);
    }
    return sum;
  }
  [[nodiscard]] T trace() const {
    if (this->column == this->row) {
      T result{0};
      for (size_t i = 0; i < this->row; i++) {
        result += this->mat[i * column + i];
      }
      return result;
    }
    throw std::runtime_error("Matrix is not square matrix");
  }
  [[nodiscard]] auto getdiagonal() {
    std::vector<T> result(this->row, 0);
    for (size_t i = 0; i < this->row; i++) {
      result[i] = this->at(i, i);
    }
    return result;
  }
  [[nodiscard]] qmatrix<T> id(size_t _row = 0) const {
    if (this->row != this->column) {
      throw std::runtime_error("Matrix is not square matrix");
    }
    if (_row == 0) {
      _row = this->row;
    }
    qmatrix<T> result(_row, _row, 0); // reverse the row and column
    for (size_t i = 0; i < _row; i++) {
      result(i, i) = 1.0;
    }
    return result;
  }
  [[nodiscard]] qmatrix<double> real() const {
    qmatrix<double> result(this->column, this->row,
                           0); // reverse the row and column
    for (size_t i = 0; i < this->dim; i++) {
      result(i) = this->at(i).real();
    }
    return result;
  }
  [[nodiscard]] qmatrix<double> imag() const {
    qmatrix<double> result(this->column, this->row,
                           0); // reverse the row and column
    for (size_t i = 0; i < this->dim; i++) {
      result(i) = this->at(i).imag();
    }
    return result;
  }
  [[nodiscard]] qmatrix<T> cTranspose() const {    // Conjugate transpose
    qmatrix<T> result(this->column, this->row, 0); // reverse the row and column
    if constexpr (std::is_same_v<T, std::complex<double>>) {
      for (size_t i = 0; i < this->row; i++) {
        for (size_t j = 0; j < this->column; j++) {
          result(j, i) = std::conj(this->mat[i * column + j]);
        }
      }
    } else {
      for (size_t i = 0; i < this->row; i++) {
        for (size_t j = 0; j < this->column; j++) {
          result(j, i) = this->mat[i * column + j];
        }
      }
    }
    return result;
  }
  [[nodiscard]] qmatrix operator*(const T &x) const {
    qmatrix result(this->row, this->column, 0);
#pragma omp parallel for // NOLINT
    for (size_t i = 0; i < this->dim; i++) {
      result(i) = this->mat.at(i) * x;
    }
    return result;
  }
  [[nodiscard]] qmatrix operator/(const T &x) const {
    qmatrix result(this->row, this->column, 0);
#pragma omp parallel for // NOLINT
    for (size_t i = 0; i < this->dim; i++) {
      result(i) = this->mat.at(i) / x;
    }
    return result;
  }
  //  void reset(const T &x = 0) {
  //    for (auto &aa : this->mat) {
  //      aa = x;
  //    }
  //  }
  T                     *data() { return this->mat.data(); }
  [[nodiscard]] const T *data() const { return this->mat.data(); }
  [[nodiscard]] auto     begin() const { return this->mat.begin(); }
  [[nodiscard]] auto     end() const { return this->mat.end(); }
  //
  //
  void                  display() {}
  [[nodiscard]] qmatrix operator+(const qmatrix<T> &rhs) const {
    // std::cout << "Started operator+";
    if (this->row == rhs.row && this->column == rhs.column) {
      qmatrix result(this->row, this->column, 0);
#pragma omp parallel for // NOLINT
      for (size_t i = 0; i < this->dim; i++) {
        result(i) = this->at(i) + rhs(i);
      }
      // std::cout << "End of operator+" << std::endl;
      return result;
    } //
      // else throw
    throw std::runtime_error("qmatrix have different size for operator+");
  }
  // Subtract
  [[nodiscard]] qmatrix operator-(const qmatrix<T> &rhs) const {
    if (this->row == rhs.row && this->column == rhs.column) {
      qmatrix result(this->row, this->column, 0);
#pragma omp parallel for // NOLINT
      for (size_t i = 0; i < this->dim; i++) {
        result(i) = this->mat[i] - rhs(i);
      }
      return result;
    }
    throw std::runtime_error("qmatrix have different size for operator -\n");
  }
  [[nodiscard]] qmatrix<T> dot(const qmatrix<T> &rhs, double talpha = 1.0) {
    if (this->column == rhs.row) {
      qmatrix result(this->row, rhs.column, 0);
      size_t  m = this->row;
      size_t  k = this->column;
      size_t  n = rhs.column;
      if (m == 0 || k == 0 || n == 0) {
        // std::cout << "One of the dimension is zero " << std::endl;
        return result;
      }
      const double alpha = talpha;
      const double beta  = 0;
      if constexpr (std::is_same_v<T, double>) {
        cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, k, alpha,
                    this->data(), k, rhs.data(), n, beta, result.data(), n);
      }
      if constexpr (std::is_same_v<T, std::complex<double>>) {
        cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, k, &alpha,
                    this->data(), k, rhs.data(), n, &beta, result.data(), n);
      }
      return result;
    }
    throw std::runtime_error("dot:qmatrix have different size for dot\n");
  }
  [[nodiscard]] std::vector<double> diag() {
    // This function diagonalizes a symmetric/harmitian matrix.
    // So eigen value are always real(double).
    // If the matrix is not symmetric the call nonsys_diag
    if (this->row != this->column) {
      throw std::runtime_error("Error: Matrix is not a square matrix! \n");
    }
    std::vector<double> w(this->row, 0);
    size_t              n = w.size();
    if (n == 0) {
      // std::cout << "This is a empty matrix" << std::endl;
      return w;
    }
    int info = -1;
    if constexpr (std::is_same_v<T, double>) {
      info = LAPACKE_dsyevd(LAPACK_ROW_MAJOR, 'V', 'U', n, this->mat.data(), n,
                            w.data());
    }
    if constexpr (std::is_same_v<T, std::complex<double>>) {
      info = LAPACKE_zheevd(
          LAPACK_ROW_MAJOR, 'V', 'U', n,
          // reinterpret_cast<__complex__ double *>(this->mat.data()), n,
          this->mat.data(), n, w.data());
    }
    // int info= LAPACKE_dsyev( LAPACK_ROW_MAJOR, 'V', 'U', n, a, n, w );
    if (info > 0) {
      std::cout << "Error:Not able to solve Eigen value problem." << std::endl;
    }
    return w;
  }
  [[nodiscard]] std::tuple<qmatrix<T>, qmatrix<T>, cm_vec>
  nonsys_diag_complex() {
    if (this->row != this->column) {
      throw std::runtime_error("Error: Matrix is not a square matrix! \n");
    }
    std::vector<T> w(this->row, 0);
    size_t         n = w.size();
    qmatrix<T>     lv(n, n, 0);
    qmatrix<T>     rv(n, n, 0);
    if constexpr (std::is_same_v<T, std::complex<double>>) {
      auto info = LAPACKE_zgeev(
          LAPACK_ROW_MAJOR, 'V', 'V', n,
          // recast the complex pointer
          // reinterpret_cast<__complex__ double *>(this->data()), n,
          this->data(), n,
          // recast the complex pointer
          w.data(),
          // recast the complex pointer
          lv.data(), n,
          // recast the complex pointer
          rv.data(), n);
/* Check for convergence */
// Normalize the vectors only for the diagonal elements
#pragma omp parallel for // NOLINT
      for (size_t i = 0; i < n; i++) {
        std::complex<double> aa{0};
        for (size_t k = 0; k < n; k++) {
          aa += std::conj(lv(k, i)) * rv(k, i);
        }
        // if (std::fabs(aa) < 1e-5) {
        //  // std::cout << "Warning:Normed:" << aa;
        //} else {
        for (size_t k = 0; k < n; k++) {
          lv(k, i) = lv(k, i) / std::conj(std::sqrt(aa));
          rv(k, i) = rv(k, i) / (std::sqrt(aa));
        }
        //}
      }
      if (info > 0) {
        throw std::runtime_error(
            "The algorithm LAPACKE_zgeev failed to compute eigenvalues.\n");
      }
    } else {
      throw std::invalid_argument(
          "nonsys_diag_complex: This function is for complex matrices");
    }
    return std::tuple(lv, rv, w);
    // return {lv, rv, w};
  }
  std::tuple<qmatrix<std::complex<T>>, qmatrix<std::complex<T>>,
             std::vector<std::complex<T>>>
  nonsys_diag_real() {
    size_t nsize = this->getrow();
    if (this->getrow() != this->getcolumn()) {
      throw std::invalid_argument(
          "nonsys_diag_real: This is not a square matrix");
    }
    if constexpr (!std::is_same_v<T, double>) {
      throw std::invalid_argument(
          "nonsys_diag_real: This is  not a qmatrix<double>! ");
    }
    std::vector<T> wr(nsize, 0);
    std::vector<T> wi(nsize, 0);
    std::vector<T> vl(nsize * nsize, 0);
    std::vector<T> vr(nsize * nsize, 0);
    //
    // timer t1("Solving exact");
    mkl_set_num_threads(200);
    auto info =
        LAPACKE_dgeev(LAPACK_ROW_MAJOR, 'V', 'V', nsize, this->data(), nsize,
                      wr.data(), wi.data(), vl.data(), nsize, vr.data(), nsize);
    // std::cout << "ExactSolver Done " << t1.getDuration() << std::endl;
    if (info > 0) {
      std::cout << "The algorithm failed to compute eigenvalues." << std::endl;
      exit(1);
    }
    std::vector<std::complex<T>> eigenvalues(nsize, 0);
    qmatrix<std::complex<T>>     leftVectors(nsize, nsize, 0);
    qmatrix<std::complex<T>>     rightVectors(nsize, nsize, 0);
// set the values
#pragma omp parallel for // NOLINT
    for (size_t j = 0; j < nsize; j++) {
      eigenvalues[j] = std::complex<T>(wr[j], wi[j]);
    }
// set eigenvector
// I dont know why the fuck this is organized this way
#pragma omp parallel for // NOLINT
    for (size_t i = 0; i < nsize; i++) {
      size_t j = 0;
      while (j < nsize) {
        if (wi[j] == static_cast<T>(0.0)) {
          leftVectors(i, j)  = vl[i * nsize + j];
          rightVectors(i, j) = vr[i * nsize + j];
          j++;
        } else {
          leftVectors(i, j) =
              std::complex<T>(vl[i * nsize + j], vl[i * nsize + j + 1]);
          leftVectors(i, j + 1) =
              std::complex<T>(vl[i * nsize + j], -vl[i * nsize + j + 1]);
          rightVectors(i, j) =
              std::complex<T>(vr[i * nsize + j], vr[i * nsize + j + 1]);
          rightVectors(i, j + 1) =
              std::complex<T>(vr[i * nsize + j], -vr[i * nsize + j + 1]);
          j += 2;
        }
      }
    }
// NOrmalize the vector
#pragma omp parallel for // NOLINT
    for (size_t i = 0; i < nsize; i++) {
      std::complex<T> aa{0};
      for (size_t k = 0; k < nsize; k++) {
        aa += std::conj(leftVectors(k, i)) * rightVectors(k, i);
      }
      if (std::fabs(aa) < 1e-5) {
        // std::cout << "Warning:Normed:" << std::fabs(aa);
      } else {
        for (size_t k = 0; k < nsize; k++) {
          leftVectors(k, i)  = leftVectors(k, i) / std::conj(std::sqrt(aa));
          rightVectors(k, i) = rightVectors(k, i) / (std::sqrt(aa));
          //  aa2 += std::conj(lv(k, i)) * rv(k, j);
        }
      }
      //       std::cout << std::endl;
    }
    return {leftVectors, rightVectors, eigenvalues};
  }
  friend std::ostream &operator<<(std::ostream &out, const qmatrix<T> &val) {
    out << "\n";
    for (size_t i = 0; i < val.row; ++i) {
      for (size_t j = 0; j < val.column; ++j) {
        out << val(i, j) << " ";
      }
      out << "\n";
    }
    return out;
  }
  qmatrix<T> krDot(const qmatrix<T> &rhs, double alpha = 1) {
    // https://en.wikipedia.org/wiki/Kronecker_product
    size_t     m = this->row;
    size_t     n = this->column;
    size_t     p = rhs.row;
    size_t     q = rhs.column;
    qmatrix<T> result(this->row * rhs.row, this->column * rhs.column, 0);
#pragma omp parallel for // NOLINT
    for (size_t r = 0; r < m; r++) {
      for (size_t s = 0; s < n; s++) {
        for (size_t v = 0; v < p; v++) {
          for (size_t w = 0; w < q; w++) {
            result((r * p + v), (s * q + w)) =
                this->at(r, s) * rhs(v, w) * alpha;
          }
        }
      }
    }
    return result;
  }
  void unitary_transform(const qmatrix<T> &eigen_vector) {
    // U^T. x . U
    // TODO(sp):
    auto result = eigen_vector.cTranspose().dot(this->dot(eigen_vector));
    *this       = result;
  }
};