Hi,

I would like to implement with Eigen (only) equivalent of expm.
Currently, I use unsupported expm feature (very fast for real) but it is very slow for complex numbers x41 slower than real where x10 or x5 could be expected.
Do you know a faster way to do with complex numbers ?

Code: Select all

 
#define EIGEN_USE_BLAS
#define EIGEN_USE_LAPACKE_STRICT

                Eigen::Map<Eigen::MatrixXd> matA((double*)A.getDataPointer(), (Eigen::Index)A.getDimensions().getRows(), (Eigen::Index)A.getDimensions().getColumns());
                Eigen::Map<Eigen::MatrixXd> matR((double*)R.getDataPointer(), (Eigen::Index)R.getDimensions().getRows(), (Eigen::Index)R.getDimensions().getColumns());
                 matR = matA.exp();
--> A = rand(500, 500);
tic(); expm(A);toc
Elapsed time is 0.232666 seconds.

                doublecomplex* Az = reinterpret_cast<doublecomplex*>((single*)A.getDataPointer());
                doublecomplex* Rz = reinterpret_cast<doublecomplex*>((single*)R.getDataPointer());
                Eigen::Map<Eigen::MatrixXcd> matA(Az, (Eigen::Index)A.getDimensions().getRows(), (Eigen::Index)A.getDimensions().getColumns());
                Eigen::Map<Eigen::MatrixXcd> matR(Rz, (Eigen::Index)R.getDimensions().getRows(), (Eigen::Index)R.getDimensions().getColumns());
                 matR = matA.exp();

--> A = rand(500, 500) + i;
tic(); expm(A);toc

Elapsed time is 9.540562 seconds.


I also tried this (but it is also slow):

Code: Select all

               Eigen::ComplexEigenSolver<Eigen::MatrixXcd> es(matA);
               auto evects = es.eigenvectors();
               auto evals = es.eigenvalues();
               for (size_t i = 0; i < static_cast<size_t>(evals.rows()); ++i)
                  evals(i) = std::exp(evals(i));
               auto evalsdiag = evals.asDiagonal();
               matR = evects * evalsdiag * evects.inverse();


Advices are welcome !!!

Thanks