0 votes
asked by (540 points)

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>
expm_small(int M, T* Hptr, T* Uptr, T tau)
  bessj = {1.266065877752008,1.13031820798497,
 bessjf = {1.,1.,1E-1,1E-2,1E-3,1E-4,1E-5,1E-6,1E-7,1E-8,1E-10,
 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;

1 Answer

0 votes
answered by (14.1k points)


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


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


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;`
Welcome to ITensor Support Q&A, where you can ask questions and receive answers from other members of the community.

Formatting Tips:
  • To format code, indent by four spaces
  • To format inline LaTeX, surround it by @@ on both sides
  • To format LaTeX on its own line, surround it by $$ above and below
  • For LaTeX, it may be necessary to backslash-escape underscore characters to obtain proper formatting. So for example writing \sum\_i to represent a sum over i.
If you cannot register due to firewall issues (e.g. you cannot see the capcha box) please email Miles Stoudenmire to ask for an account.

To report ITensor bugs, please use the issue tracker.