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