How to use Matrix class

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;
}
}


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

commented by (540 points)
Thanks Matt, i'll take a look and try.
commented by (540 points)
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;
commented by (70.1k points)
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.
commented by (540 points)
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;