Fast matrix exponential of a complex symmetric tridiagonal matrix

Question

Basically I need the above. I have trawled the Googles and cannot find a way of implementing this.

I've found this function here http://www.guwi17.de/ublas/examples/ but it is too slow. I even wrote my own Pade Approximation following MATLAB's routine but it was only fractionally faster than the one in the link.

What stuns me is how quick Mathematica is able to compute matrix exponentials (whether it cares if the matrix is tridiag I don't know).

Is anyone able to help?

EDIT: Here's what I've come up with, any comments? Will hopefully be useful for future readers

I've been out of c++ for a while so the below code may be a bit scrappy / slow, so please enlighten me if you see improvements.

//Program will compute the matrix exponential expm( i gA ). where i = sqrt(-1) and gA is a matrix

#include <iostream>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>


int main() {

    int n = 100;

    unsigned i = 0;
    unsigned j = 0;

    gsl_complex img = gsl_complex_rect (0.,1.);

    gsl_matrix * gA = gsl_matrix_alloc (n, n);


//Set gA:       
    for ( i = 0; i<n; i++) {
        for ( j = 0; j<=i; j++) {
            if( i == j ) {
                if( i == 0 || i == n-1 ) {
                    gsl_matrix_set (gA, i, j, 1.);
                } else {
                    gsl_matrix_set (gA, i, j, 2.);
                }
            } else if( j == i-1 ) {
                gsl_matrix_set (gA, i, j, 1.);
                gsl_matrix_set (gA, j, i, 1.);
            } else {
                gsl_matrix_set (gA, i, j, 0.);
                gsl_matrix_set (gA, j, i, 0.);
        }
    }
}


//get SVD of gA 
    gsl_matrix * gA_t = gsl_matrix_alloc (n, n);
    gsl_matrix_transpose_memcpy (gA_t, gA);                 

    gsl_matrix * U = gsl_matrix_alloc (n, n);
    gsl_matrix * V= gsl_matrix_alloc (n, n);
    gsl_vector * S = gsl_vector_alloc (n);


    gsl_vector * work = gsl_vector_alloc (n);
    gsl_linalg_SV_decomp (gA_t, V, S, work);
    gsl_vector_free(work);

    gsl_matrix_memcpy (U, gA_t);


//Take exponential of S and form matrix
    gsl_matrix_complex * Sp = gsl_matrix_complex_alloc (n, n);
    gsl_matrix_complex_set_zero (Sp);
    for (i = 0; i < n; i++) {
        gsl_matrix_complex_set (Sp, i, i, gsl_complex_exp ( gsl_complex_mul_real ( img, gsl_vector_get(S, i) ) ) ); // Vector 'S' to matrix 'Sp'
    }

    gsl_matrix_complex * Uc = gsl_matrix_complex_alloc (n, n);


//convert U to a complex matrix for next step   
    for (i = 0; i < n; i++) {
        for ( j = 0; j<n; j++) {
            gsl_matrix_complex_set (Uc, i, j, gsl_complex_rect ( gsl_matrix_get(U, i, j), 0. ) );       
        }
    }


//recombine U S Utranspose  
    gsl_matrix_complex * tA = gsl_matrix_complex_alloc (n, n);
    gsl_matrix_complex * eA = gsl_matrix_complex_alloc (n, n);
    gsl_blas_zgemm (CblasNoTrans, CblasNoTrans, GSL_COMPLEX_ONE, Uc, Sp, GSL_COMPLEX_ZERO, tA);
    gsl_blas_zgemm (CblasNoTrans, CblasTrans, GSL_COMPLEX_ONE, tA, Uc, GSL_COMPLEX_ZERO, eA);


//Print result  
    std::cout << "eA:" << std::endl;
    for (i = 0; i < n; i++)  
        for (j = 0; j < n; j++)
            printf ("%g + %g i\n", GSL_REAL(gsl_matrix_complex_get (eA, i, j)), GSL_IMAG(gsl_matrix_complex_get (eA, i, j)));
    std::cout << "\n" << std::endl;


//Free up
    gsl_matrix_free(gA_t);
    gsl_matrix_free(U);
    gsl_matrix_free(gA);
    gsl_matrix_free(V);
    gsl_vector_free(S);
    gsl_matrix_complex_free(Uc);
    gsl_matrix_complex_free(tA);    
    gsl_matrix_complex_free(eA);    

    return 0;
}        

Answer 1

The code above I've posted into my question works very well. Another improvement I found was to scale the columns in a copy of V by the eigenvalues rather than perform a full matrix multiplication. With zgemm being the slowest part of this code, removing one of the zgemm functions sped it up nearly two-fold. I'll post a full code listing here shortly.