Complex matrix exponential in C++

Question

Is it actually possible to calculate the Matrix Exponential of a Complex Matrix in c / c++?

I've managed to take the product of two complex matrices using blas functions from the GNU Science Library. for matC = matA * matB:

gsl_blas_zgemm (CblasNoTrans, CblasNoTrans, GSL_COMPLEX_ONE, matA, matB, GSL_COMPLEX_ZERO, matC);

And I've managed to get the matrix exponential of a matrix by using the undocumented

gsl_linalg_exponential_ss(&m.matrix, &em.matrix, .01);

But this doesn't seems to accept complex arguments.

Is there anyway to do this? I used to think c++ was capable of anything. Now I think its outdated and cryptic...

See example at: http://www.guwi17.de/ublas/examples/ (the bottom of the page). — Jerry Coffin, Apr 06 '12 at 21:01
This problem has nothing to do with the C++ language itself, but with a GSL library, so don't dislike C++ for that :) Also, there are no tools capable of anything. — Rafał Rawicki, Apr 06 '12 at 21:02
Thanks Jerry, I'm now using the function given at the bottom of that page. It works, it's unfortunately a little slower than I would have like but it looks like I may just be able to take fewer sample points and run some interpolation (which is much quicker). Many thanks for pointing that out! — Quantum_Oli, Apr 10 '12 at 14:12

score 3 · Answer 1 · answered Apr 12 '12 at 22:19

Several options:

modify the gsl_linalg_exponential_ss code to accept complex matrices
write your complex NxN matrix as real 2N x 2N matrix
Diagonalize the matrix, take the exponential of the eigenvalues, and rotate the matrix back to the original basis
Using the complex matrix product that is available, implement the matrix exponential according to it's definition: exp(A) = sum_(n=0)^(n=infinity) A^n/(n!)

You have to check which methods are appropriate for your problem.

C++ is a general purpose language. As mentioned above, if you need specific functionality you have to find a library that can do it or implement it yourself. Alternatively you could use software like MatLab and Mathematica. If that's too expensive there are open source alternatives, e.g. Sage and Octave.

(2) wastes some memory, but is probably your best option, since (3) will only work if you know that your matrix is normal (AA^* = A^*A), (4) is not the best way to implement exponential (use exp(A / 2^n)^(2^n) = exp (A), find n such that the norm of A / 2^n is small and taylor expand exp (A/2^n) = B. Now compute B^(2^n) with n multiplications) — Alexandre C., Feb 04 '13 at 20:32

score 1 · Answer 2 · answered Apr 06 '12 at 21:04

"I used to think c++ was capable of anything" - if a general-purpose language has built-in complex math in its core, then something is wrong with that language.

Fur such very specific tasks there is a well-accepted solution: libraries. Either write your own, or much better, use an already existing one.

I myself rarely need complex matrices in C++, I always used Matlab and similar tools for that. However, this http://www.mathtools.net/C_C__/Mathematics/index.html might be of interest to you if you know Matlab.

There are a couple other libraries which might be of help:

score 0 · Answer 3 · edited Feb 04 '13 at 18:11

I was also thinking to do the same, writing your complex NxN matrix as real 2N x 2N matrix is the best way to solve the problem, then use gsl_linalg_exponential_ss().

Suppose A=Ar+i*Ai, where A is the complex matrix and Ar and Ai are the real matrices. Then write the new matrix B=[Ar Ai ;-Ai Ar] (Here the matrix is written in matlab notation). Now calculate the exponential of B, that is eB=[eB1 eB2 ;eB3 eB4], then exponential of A is given by, eA=eB1+1i.*eB2
(summing the matrices eB1 and 1i.*eB2).

score 0 · Answer 4 · edited Feb 05 '13 at 23:13

I have written a code to calculate the matrix exponential of the complex matrices with the gsl function, gsl_linalg_exponential_ss(&m.matrix, &em.matrix, .01); Here you have the complete code, and the compilation results. I have checked the result with the Matlab and result agrees.

#include <stdio.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
void my_gsl_complex_matrix_exponential(gsl_matrix_complex *eA, gsl_matrix_complex *A, int dimx)
{
    int j,k=0;
    gsl_complex temp;
    gsl_matrix *matreal =gsl_matrix_alloc(2*dimx,2*dimx);
    gsl_matrix *expmatreal =gsl_matrix_alloc(2*dimx,2*dimx);
    //Converting the complex matrix into real one using A=[Areal, Aimag;-Aimag,Areal]
    for (j = 0; j < dimx;j++)
        for (k = 0; k < dimx;k++)
        {
            temp=gsl_matrix_complex_get(A,j,k);
            gsl_matrix_set(matreal,j,k,GSL_REAL(temp));
            gsl_matrix_set(matreal,dimx+j,dimx+k,GSL_REAL(temp));
            gsl_matrix_set(matreal,j,dimx+k,GSL_IMAG(temp));
            gsl_matrix_set(matreal,dimx+j,k,-GSL_IMAG(temp));
        }

    gsl_linalg_exponential_ss(matreal,expmatreal,.01);

    double realp;
    double imagp;
    for (j = 0; j < dimx;j++)
        for (k = 0; k < dimx;k++)
        {
            realp=gsl_matrix_get(expmatreal,j,k);
            imagp=gsl_matrix_get(expmatreal,j,dimx+k);
            gsl_matrix_complex_set(eA,j,k,gsl_complex_rect(realp,imagp));
        }
    gsl_matrix_free(matreal);
    gsl_matrix_free(expmatreal);
}

int main()
{
    int dimx=4;
    int i, j ;
    gsl_matrix_complex *A = gsl_matrix_complex_alloc (dimx, dimx);
    gsl_matrix_complex *eA = gsl_matrix_complex_alloc (dimx, dimx);

    for (i = 0; i < dimx;i++)
    {
        for (j = 0; j < dimx;j++)
        {
            gsl_matrix_complex_set(A,i,j,gsl_complex_rect(i+j,i-j));
            if ((i-j)>=0)
            printf("%d+%di ",i+j,i-j);
            else
            printf("%d%di  ",i+j,i-j);

        }
        printf(";\n");
    }
    my_gsl_complex_matrix_exponential(eA,A,dimx);
    printf("\n Printing the complex matrix exponential\n");
    gsl_complex compnum;
    for (i = 0; i < dimx;i++)
    {
        for (j = 0; j < dimx;j++)
        {
            compnum=gsl_matrix_complex_get(eA,i,j);
            if (GSL_IMAG(compnum)>=0)
                printf("%f+%fi\t ",GSL_REAL(compnum),GSL_IMAG(compnum));
            else
                printf("%f%fi\t ",GSL_REAL(compnum),GSL_IMAG(compnum));
        }
        printf("\n");
    }
    return(0);
}

Complex matrix exponential in C++

4 Answers4

Linked