Algebra.cpp

#include <iostream>

#include "Algebra.hpp"

namespace Algebra{

  using namespace Eigen;

  QMat::QMat(const QuMatrix& m){
    mat = m;
  }

  QMat::QMat(int r, int c){
    QuMatrix m(r,c);

    for(int i=0; i<r; i++){
      for(int j=0; j<c; j++){
        m(i,j) = 0;
      }
    }

    mat = m;
  }

  QMat::QMat(int n, bool ident /*= false*/){
    QuMatrix m(n,n);

    for(int i=0; i<n; i++){
      for(int j=0; j<n; j++){
        m(i,j) = 0;
      }
    }

    if(ident == true){
      for(int o=0; o<n; o++){
        m(o,o) = 1; 
      }
    }

    mat = m;
  }

  QMat::QMat(std::vector<std::vector<std::complex<double>>> l){
    QuMatrix m(l.size(), l[0].size()); 

    for(int i=0; i<l.size(); i++){
      for(int j=0; i<l[0].size(); i++){
        m(i,j) = l[i][j];
      }
    }

    mat = m;
  }

  QMat QMat::operator+(const QMat& m){
    return QMat(this->mat + m.mat);
  }

  QMat QMat::operator-(const QMat& m){
    return QMat(this->mat - m.mat);
  }

  QMat QMat::operator*(const QMat& m){
    return QMat(this->mat * m.mat);
  }

  QMat QMat::operator*(const std::complex<double> c){
    return QMat(this->mat * c);
  }

  void QMat::operator+=(const QMat& m){
    this->mat+=m.mat;
  }

  void QMat::operator-=(const QMat& m){
    this->mat-=m.mat;
  }

  void QMat::operator*=(const QMat& m){
    this->mat*=m.mat;
  }

  void QMat::operator*=(const std::complex<double> c){
    this->mat*=c;
  }

  void QMat::setElement(int i, int j, std::complex<double> c){
    this->mat(i,j) = c; 
  }

  std::complex<double> QMat::getElement(int i, int j){
    return this->mat(i,j);
  }

  int QMat::rows(){
    return this->mat.rows();
  }

  int QMat::cols(){
    return this->mat.cols();
  }

  std::complex<double> QMat::trace(){
    return this->mat.trace();
  }

  std::complex<double> QMat::inner_product(const QMat& m){
    assert(this->mat.cols() == 1);
    assert(m.mat.cols() == 1);
    assert(this->mat.rows() == m.mat.rows());
    std::complex<double> sum = 0;
    for(int i=0; i<this->mat.rows(); i++){
      sum += this->mat(i,0) * std::conj(m.mat(i,0));
    }
    return sum;
  }

  QMat QMat::conjugateTranspose(){
    return QMat(mat.adjoint());
  }

  QMat QMat::orthoNormalize(){
    return QMat(mat.householderQr().householderQ());
  }

  QMat QMat::tensor(const QMat& m){
    return QMat(kroneckerProduct(this->mat, m.mat));
  }

  QMat QMat::eigenValues(){
    ComplexEigenSolver<QuMatrix> ces(mat, false);
    return QMat(ces.eigenvalues());
  }

  QMat QMat::eigenVectors(){
    ComplexEigenSolver<QuMatrix> ces(mat);
    return QMat(ces.eigenvectors());
  }

  void QMat::serialize(){
    for(int i=0; i < this->mat.rows(); i++){
      for(int j=0; j < this->mat.cols(); j++){
        if(j != this->mat.cols() - 1){
          std::cout << this->mat(i,j) << ",";
        }else{
          std::cout << this->mat(i,j);
        }
      }
      std::cout << ";" << std::endl;
    }
  }

  std::ostream& operator<<(std::ostream& out, const QMat& m){
    return out << m.mat;
  }
};