Boost logo

Ublas :

From: Tsai Dung-Bang (dbtsai_at_[hidden])
Date: 2007-11-22 10:12:54

Dear all

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
ublas library.

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
real 3m49.192s
user 3m48.926s
sys 0m0.264s,
and in Fortran code, it only took
real 0m16.886s
user 0m16.885s
sys 0m0.000s.

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


Best Wish

Tsai, Dung-Bang
National Taiwan University, Department of physics