How to use Matrix class

Question

hi,
I'm trying to implement Chebyshev expansion matrix exponential (see P.6 @paper) under itensor but lack of documents about Matrix class. The algorithm realized by Eigen is as follows and is already working, can any one implement this code under pure itenor framework. The matrix and result are both stored in std::vector<Real/Cplx>

std::vector<Real/Cplx> H stores matrix H in column major, and the result e^(tau*H) is then returned in std::vector<Real/Cplx> U, M is the dimension of matrix H.
usage Eigen::expm_small(M,H.data(),U.data(),tau);

//Chebyshev expansion ExpM U = e^(tau*H)
template<typename T>
void
expm_small(int M, T* Hptr, T* Uptr, T tau)
{
  std::vector<double>
  bessj = {1.266065877752008,1.13031820798497,
         2.714953395340766,4.43368498486638,
         5.474240442093732,5.429263119139439,
         4.497732295429515,3.19843646240199,
         1.992124806672796,1.103677172551734,
         5.505896079673748,2.497956616984982,
         1.03915223067857,3.991263356414401,
         1.423758010825657,4.740926102561494,
         1.480180057208297,4.349919494944169,
         1.207428927279753,3.175356737059445};
 std::vector<double>
 bessjf = {1.,1.,1E-1,1E-2,1E-3,1E-4,1E-5,1E-6,1E-7,1E-8,1E-10,
          1E-11,1E-12,1E-14,1E-15,1E-17,1E-18,1E-20,1E-21,1E-23};
 Map<Matrix<T,Dynamic,Dynamic>> hmat(Hptr, M, M);
 Map<Matrix<T,Dynamic,Dynamic>> umat(Uptr, M, M);
 Matrix<T,Dynamic,Dynamic> t0 = Matrix<T,Dynamic,Dynamic>::Identity(M,M);
 Matrix<T,Dynamic,Dynamic> t1 = hmat * tau;
 Matrix<T,Dynamic,Dynamic> m = t1;
 umat = t0 * bessj[0];
 umat += t1 * bessj[1];
 for(int i = 2; i < 19; ++i)
 {
   Matrix<T,Dynamic,Dynamic> tk = 2 * (m * t1) - t0;
   umat +=  (tk * bessj[i]) * bessjf[i];
   t0 = t1;
   t1 = tk;
 }
}

Answer 1

Hello,

You can take a look at the unit tests for the ITensor Matrix class to get an idea for how it is used:

https://github.com/ITensor/ITensor/blob/v3/unittest/matrix_test.cc#L922

Let us know if you have questions/some functionality that you need is not there.

Cheers,
Matt

Comments

Thanks Matt, i'll take a look and try.

---

Hi, Matt
Basically the code works, but only for Real tau. Will you please look into this code (https://gist.github.com/empter/0ffe892d7b3316c6f47b47b094592d16) and figure out what's the problem with Cplx tau.
The error is `error: no match for ‘operator*’ (operand types are ‘itensor::TenRef<itensor::MatRangeT<0>, std::complex<double> >’ and ‘std::complex<double>’)
   Mat<T> t1 = subMatrix(H,0,M,0,M) * tau;`

---

Hi, could you try first saving the subMatrix to a Mat<T> t1, then multiplying t1 in-place by tau?

There are some limitations to the Matrix layer of ITensor, as I never completely "polished" it for use by end users, but only designed it as an internal developer-level thing since none of the other matrix libraries in C++ fit our purposes.

---

Hi, Miles
I tried `Mat<T> t1 = subMatrix(H,0,M,0,M);`
still error `error: conversion from ‘itensor::TenRef<itensor::MatRangeT<0>, std::complex<double> >’ to non-scalar type ‘itensor::Mat<std::complex<double> > {aka itensor::Ten<itensor::MatRangeT<0>, std::complex<double> >}’ requested
   Mat<T> t1 = subMatrix(H,0,M,0,M);`
is there any efficient way to save the subMatrix to a Mat<T> t1?

And the matrix class seems not supporting *= operator for complex scalars, using `H *= tau;` directly gives error `error: no match for ‘operator*=’ (operand types are ‘itensor::Mat<std::complex<double> > {aka itensor::Ten<itensor::MatRangeT<0>, std::complex<double> >}’ and ‘std::complex<double>’)
   H *= tau;`