From: Tsai Dung-Bang (dbtsai_at_[hidden])
Date: 2007-11-22 10:12:54
As far as everyone know, matrix exponential is very important in
physics and engineering.
For example, a differential equation , d/dt v = M v where v is vector
and M is matrix could be
solved by v(t) = v(0)*exp(Mt), where v(0) is the initial state. Such
kind of equations are
often disappeared in lots of problem. Also, in quantum theory,
Schrodinger equation has the same form, ihd/dt V = H V, where V is the
state and H is Hamiltonian, and in most case, we could use matrix
representation such that H could be written as matrix, and state
could be written as vector.
As a result it's worthwhile to include this function in the mainline
In my library, see attachment, I implement it using Pade approximation
where matlab also uses the same algorithm. In my quantum system, it
works great, and I have compare the result between matlab and Fortran
version with random matrix, they have the same result.
But the performance is not as good as Fortran version, which is name
after expokit package. I use the same algorithm as expokit one, but
with the same 1e6 times 5 by 5 random matrix operation, in my C++
code, it took
and in Fortran code, it only took
I have no idea about that, I wish someone could investigate it.
Also, I have some question about template programming, the following
code is took from my code.
template<typename MATRIX> MATRIX expm_pad(const MATRIX &H, const int p = 6)
typedef typename MATRIX::value_type value_type;
typedef typename MATRIX::size_type size_type;
typedef double real_value_type; // Correct me. Need to modify.
assert(H.size1() == H.size2());
In the typedef real_value_type, how could I do this?
if( value_type is complex)
typedef typename value_type::value_type real_value_type;
typedef value_type real_value_type
National Taiwan University, Department of physics