[ublas] Matrix exponential library using Pade approximation

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;
else
  typedef value_type real_value_type

Thanks

Best Wish

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