|
|
Ublas : |
From: Riccardo Rossi (rrossi_at_[hidden])
Date: 2007-11-22 10:33:13
Hi Tsai,
try to use a "bounded_matrix" to my knowledge it is faster than the
"matrix" (and it should be given that everything is known at compile time!!)
good luck
Riccardo
Tsai Dung-Bang wrote:
> 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
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> ublas mailing list
> ublas_at_[hidden]
> http://lists.boost.org/mailman/listinfo.cgi/ublas