Program Listing for File matrix.cpp

↰ Return to documentation for file (src_util/matrix.cpp)

#include <Eigen/Dense>

#include "matrix.hpp"

using namespace std;

/*
 NOTE : The general convention of this code is to use "Column Major Matrices",
 stored in linear arrays, such that if A is a r x c matrix, then

 A[i][j] --->  A[i+r*j]

 Eigen also uses column-major storage by default, so Eigen::Map wraps
 the v.data() buffer directly with zero copy.
 */

template<> void matrix<double>::inverse()
{
    QCM_ASSERT(r==c);
    if(r==0) return;
    if(r==1){v[0] = 1/v[0]; return;}
    Eigen::Map<Eigen::MatrixXd> M(v.data(), r, c);
    M = M.inverse().eval();
}


template<> bool matrix<Complex>::is_real(double accuracy)
{
    for(size_t i=0; i<v.size(); ++i) if(abs(v[i].imag()) > accuracy) return false;
    return true;
}


template<> void matrix<double>::eigensystem(vector<double> &d, matrix<double> &U) const
{
    QCM_ASSERT(r == c);
    QCM_ASSERT(U.r == r and U.c == c);
    QCM_ASSERT(d.size() == r);

    Eigen::Map<const Eigen::MatrixXd> A(v.data(), r, c);
    Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> solver(A);
    if(solver.info() != Eigen::Success)
    {
        cout << "dimension : " << r << " x " << c << endl;
        cout << *this << endl;
        throw std::runtime_error("SelfAdjointEigenSolver failed (eigensystem<double>)");
    }
    Eigen::Map<Eigen::VectorXd>(d.data(), r) = solver.eigenvalues();
    Eigen::Map<Eigen::MatrixXd>(U.v.data(), r, c) = solver.eigenvectors();
}


template<> void matrix<double>::eigenvalues(vector<double> &d)
{
    QCM_ASSERT(r == c);
    Eigen::Map<const Eigen::MatrixXd> A(v.data(), r, c);
    Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> solver(A, Eigen::EigenvaluesOnly);
    QCM_ASSERT(solver.info() == Eigen::Success);
    Eigen::Map<Eigen::VectorXd>(d.data(), r) = solver.eigenvalues();
}


template<> void matrix<Complex>::inverse()
{
    QCM_ASSERT(r==c);
    if(r==0) return;
    Eigen::Map<Eigen::MatrixXcd> M(v.data(), r, c);
    M = M.inverse().eval();
}


template<> Complex matrix<Complex>::determinant()
{
    if(r==0) return Complex(0);
    Eigen::Map<Eigen::MatrixXcd> M(v.data(), r, c);
    return M.determinant();
}


template<> void matrix<Complex>::eigensystem(vector<double> &d, matrix<Complex> &U) const
{
    QCM_ASSERT(r == c);

    Eigen::Map<const Eigen::MatrixXcd> A(v.data(), r, c);
    Eigen::SelfAdjointEigenSolver<Eigen::MatrixXcd> solver(A);
    if(solver.info() != Eigen::Success)
    {
        cout << "dimension : " << r << " x " << c << endl;
        cout << *this << endl;
        throw std::runtime_error("SelfAdjointEigenSolver failed (eigensystem<complex>)");
    }
    Eigen::Map<Eigen::VectorXd>(d.data(), r) = solver.eigenvalues();
    Eigen::Map<Eigen::MatrixXcd>(U.v.data(), r, c) = solver.eigenvectors();
}


template<> void matrix<Complex>::eigenvalues(vector<double> &d)
{
    QCM_ASSERT(r == c);
    Eigen::Map<const Eigen::MatrixXcd> A(v.data(), r, c);
    Eigen::SelfAdjointEigenSolver<Eigen::MatrixXcd> solver(A, Eigen::EigenvaluesOnly);
    QCM_ASSERT(solver.info() == Eigen::Success);
    Eigen::Map<Eigen::VectorXd>(d.data(), r) = solver.eigenvalues();
}


template<> double matrix<double>::norm()
{
    Eigen::Map<const Eigen::MatrixXd> M(v.data(), r, c);
    return M.norm();
}


template<> double matrix<Complex>::norm()
{
    Eigen::Map<const Eigen::MatrixXcd> M(v.data(), r, c);
    return M.norm();
}


template<> double matrix<double>::norm2()
{
    Eigen::Map<const Eigen::MatrixXd> M(v.data(), r, c);
    return M.squaredNorm();
}


template<> double matrix<Complex>::norm2()
{
    Eigen::Map<const Eigen::MatrixXcd> M(v.data(), r, c);
    return M.squaredNorm();
}


template<> bool matrix<double>::is_orthogonal(double accuracy)
{
    if(r!=c) return false;
    Eigen::Map<const Eigen::MatrixXd> M(v.data(), r, c);
    Eigen::MatrixXd diff = M.transpose() * M - Eigen::MatrixXd::Identity(r, r);
    return diff.cwiseAbs().maxCoeff() <= accuracy;
}


template<> bool matrix<double>::cholesky(matrix<double> &A)
{
    Eigen::Map<const Eigen::MatrixXd> Ma(A.v.data(), A.r, A.c);
    Eigen::LLT<Eigen::MatrixXd> llt(Ma);
    if(llt.info() != Eigen::Success) return false;
    set_size(A.r);
    Eigen::Map<Eigen::MatrixXd>(v.data(), r, c) = llt.matrixL();
    return true;
}


template<> bool matrix<double>::triangular_inverse()
{
    Eigen::Map<Eigen::MatrixXd> M(v.data(), r, c);
    Eigen::MatrixXd Linv = M.triangularView<Eigen::Lower>().solve(
        Eigen::MatrixXd::Identity(r, r));
    M = Linv;
    return true;
}



/*
transforms a real matrix into a complex one
**/
matrix<Complex> to_complex_matrix(const matrix<double> &x) {
    matrix<Complex> xc(x.r, x.c);
    for (size_t i = 0; i < x.v.size(); i++) xc.v[i] = x.v[i];
    return xc;
}

/*
Builds a hermitian complex matrix from a real vector
**/
matrix<Complex> hermitian_matrix_from_real_vector(size_t d, const vector<double> &x){
    matrix<Complex> M(d);
    size_t k = 0;
    for(int i = 0; i<d; i++){
      M(i,i) = x[k++];
      for(int j = 0; j<i; j++){
        M(i,j) = Complex(x[k], x[k+1]);
        k += 2;
        M(j,i) = conjugate(M(i,j));
      }
    }
    return M;
}

/*
Populates a real vector from a hermitian complex matrix
**/
void hermitian_matrix_to_real_vector(const matrix<Complex> &M, double *x){
    size_t k = 0;
    for(int i = 0; i<M.r; i++){
      x[k++] = M(i,i).real();
      for(int j = 0; j<i; j++){
        x[k++] = M(i,j).real();
        x[k++] = M(i,j).imag();
      }
    }
}