.. _program_listing_file_src_util_matrix.cpp:

Program Listing for File matrix.cpp
===================================

|exhale_lsh| :ref:`Return to documentation for file <file_src_util_matrix.cpp>` (``src_util/matrix.cpp``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: 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();
         }
       }
   }