[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: Using expm in C++

From: loisp
Subject: Re: Using expm in C++
Date: Sat, 11 Feb 2012 19:02:42 -0800 (PST)

Hi Matyas,

First, thanks! That's a lot of effort. 

I haven't gone through your entire reply carefully just yet but I thought I should tell you what I did. I wrote a Standalone C++ program, used Octave's eig function to do the job (with the help of this document) and compiled with mkoctfile. I'm pasting the relevant snippet below.

EIG eig = EIG(hamiltonian);                                /*Eigensystem of Hamiltonian*/

ComplexColumnVector eigen = eig.eigenvalues();

D = ComplexMatrix(ComplexDiagMatrix(eigen));

ComplexMatrix U = (eig.eigenvectors());     /*Eigenvectors*/

ComplexMatrix V = U.inverse();              /*Eigenvectors_inverse*/

for(int i = 0; i < N+1; i++) 


   for(int j = 0; j < N+1; j++) 




       D (i, j) = exp(z*D(i,j));




time_evolution = ComplexMatrix (N+1, N+1);  /*Time evolution operator U*/

time_evolution = U*D*V;

I checked my results and it seems to work pretty well. Since I don't use MATLAB I steered clear of mex-files. 

The matrix I am using is a representation of the Hamiltonian operator and it's Hermitian (it is not tridiagonal - the diagonal elements are zero). I think it satisfies Wikipedia's definition of a Hamiltonian matrix - it's symmetric (real counterpart of Hermitian), with real eigen values and trace zero. Please let me know if I'm missing something. 
My (computational bit) problem was getting the time-evolution matrix (which is an exponential function of the Hamiltonian). 


On Sun, Feb 12, 2012 at 6:03 AM, Matyas Sustik [via Octave] <[hidden email]> wrote:
Hi Lois,

I found your first note on this:

I'm writing a C++ program and want to use the Octave function expm()
(exponential of a matrix). Here is the relevant snippet from my code.

       hamiltonian = Matrix (N+1, N+1);
       for(int i = 0; i < N+1; i++)
               for(int j = 0; j < N+1; j++)
                       if((i+1 == j) || (i == j+1))
                               hamiltonian (i, j) = -1;

I commented earlier that the eigendecomposition could be used.  Below is the,
code I came up with:

// hamexpm.C

#include <math.h>
#include <stddef.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdint.h>
#ifdef GDEBUG
#include "startgdb.h"

// It would be preferable to use an include such as lapack.h.  Except
// lapack.h is not available from the octave or liblapack-dev
// packages.  You may find one in the extern/include directory of your
// Matlab installation.
extern "C" {
    // Matrix multiplication:
    int dgemm_(const char*, const char*, ptrdiff_t*, ptrdiff_t*,
           ptrdiff_t*, double*, double*, ptrdiff_t*, double*,
           ptrdiff_t*, double*, double*, ptrdiff_t*);
    // Eigendecomposition of a tridiagonal matrix using the state of
    // the art RRR method.
    void dstemr_(char* jobz,
         char* range,
         ptrdiff_t* n,
         double* d,
         double* e,
         double* vl,
         double* vu,
         ptrdiff_t* il,
         ptrdiff_t* iu,
         ptrdiff_t* m,
         double* w,
         double* z,
         ptrdiff_t* ldz,
         ptrdiff_t* nzc,
         ptrdiff_t* isuppz,
         ptrdiff_t* tryrac,
         double* work,
         ptrdiff_t* lwork,
         ptrdiff_t* iwork,
         ptrdiff_t* liwork,
         ptrdiff_t* info);

extern "C" void hamexpm(ptrdiff_t& n, double* A)
#ifdef GDEBUG

    double* d = (double*) malloc(n*sizeof(double));
    memset(d, 0, n*sizeof(double));
    double* e = (double*) malloc(n*sizeof(double));
    for (ptrdiff_t i = 0; i < n; i++)
    e[i] = -1.0;
    ptrdiff_t m = 0;
    double* w = (double*) malloc(n*sizeof(double));
    double* z = (double*) malloc(n*n*sizeof(double));
    ptrdiff_t* isuppz = (ptrdiff_t*) malloc(2*n*sizeof(double));
    ptrdiff_t tryrac = 1;
    double* work = (double*) malloc(18*n*sizeof(double));
    ptrdiff_t lwork = 18*n;
    ptrdiff_t* iwork = (ptrdiff_t*) malloc(10*n*sizeof(double));
    ptrdiff_t liwork = 10*n;
    ptrdiff_t info = 0;
//    dstemr_((char*) "V", (char*) "A", &n, d, e, 0, 0, 0, 0, &m, w, z, &n,
//        &n, &tryrac, work, &lwork, iwork, &liwork, &info); 

    dstemr_((char*) "V", (char*) "A", &n, d, e, 0, 0, 0, 0, &m, w, z, &n,
        &n, isuppz, &tryrac, work, &lwork, iwork, &liwork, &info);

    if (info != 0)
    printf("Ooops, dstemr() failed.  That was not supposed "
           "to happen...\n");

    // I like to free memory as sson as possible...
    // A = z*exp(w)*z'; Using BLAS the following could be faster.
    // Symmetry could also be exploited, etc.
    double* z1 = (double*) malloc(n*n*sizeof(double));
    memcpy(z1, z, n*n*sizeof(double));
    for (ptrdiff_t i = 0; i < n; i++) {
    double x = exp(w[i]);
    for (ptrdiff_t j = 0; j < n; j++)
        z1[i*n + j] = z[i*n + j]*x;

    double alpha = 1.0;
    double beta = 0.0;
    dgemm_("N", "T", (ptrdiff_t*)&n, (ptrdiff_t*)&n,
       (ptrdiff_t*)&n, &alpha, z1, (ptrdiff_t*)&n, z,
       (ptrdiff_t*)&n, &beta, A, (ptrdiff_t*)&n);


I also have a wrapper file I used to create a MEX object (shame on me I still have not learned how to make OCT files):

// hamexpm-mex.C

// This is the MEX wrapper for hemexpm.  The algorithm is in hemexpm.C.

// Invocation from within Matlab or Octave:
// [A] = hamexpm(n)

// A = expm(H), where H is nxn with A(i,i+1) = A(i+1,i) = -1,
// other entries are 0-s.

#include <mex.h>
#include <stddef.h>
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <stdint.h>

extern "C" void hamexpm(ptrdiff_t& n, double* A);

#define EPS (double(2.22E-16))

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
    ptrdiff_t n = mxGetScalar(prhs[0]);
    plhs[0] = mxCreateDoubleMatrix(n, n, mxREAL);
    hamexpm(n, mxGetPr(plhs[0]));

And for your convenience the Makefile I used:

# Makefile

OCTAVE_INCLUDES = -I/usr/include/octave
# Octave has a nice feature to tell us the library locations:
OCTAVE_LIBS = $(shell mkoctfile -p LFLAGS)


CXXFLAGS = -Wall -fpic -pthread -shared  -fno-omit-frame-pointer -ansi -D_GNU_SOURCE -D_FILE_OFFSET_BITS=64

LDOCTMEXFLAGS = -shared -Wl,-Bsymbolic -loctave -lcruft -llapack -lblas -lfftw3 -lfftw3f -lreadline -lncurses -ldl -lhdf5 -lz -lm -lgfortranbegin -lgfortran $(OCTAVE_LIBS)


# C++ rules for MEX wrappers:
%-oct_g.o: %-mex.C
    $(CXX) $(CXXDBGFLAGS) $(OCTAVE_INCLUDES) -c $< -o $@

%-oct.o: %-mex.C
    $(CXX) $(CXXOPTFLAGS) $(OCTAVE_INCLUDES) -c $< -o $@

# C++ rules for algorithm programs:
%_g.o: %.C
    $(CXX) $(CXXDBGFLAGS) -c $< -o $@

%.o: %.C
    $(CXX) $(CXXOPTFLAGS) -c $< -o $@

# Link Octave executables:
%_g.mex : %-oct_g.o %_g.o
    g++ $(LDOCTMEXFLAGS) $^ -o $@

%.mex : %-oct.o %.o
    g++ $(LDOCTMEXFLAGS) $^ -o $@


clean :
    rm -f *.o *.oct *.mex

So just issue make hamexpm.mex to create the executable.  I use linux and not OSX, I hope you can get this to work.

Please let me know if this was helpful. I would also be interested in hearing some more details about this particular problem.  On Wikipedia Hamiltonian matrices are defined differently than your use.  Can you point
me to any details on why you call the tridiagonal matrix a Hamiltonian?

Take care,

P.S.: You may remove the reference to startgdb(), I do not want to edit the above... You do not need that for the optimized version, and ifdefs make sure it does not interfere with optimized compile.  I use a little utility that pops up a debugger if you compile and run hamexp_g.mex.  I guess others debug C code under Octave using some similar trick. M.

Help-octave mailing list
[hidden email]

If you reply to this email, your message will be added to the discussion below:
To unsubscribe from Using expm in C++, click here.

View this message in context: Re: Using expm in C++
Sent from the Octave - General mailing list archive at

reply via email to

[Prev in Thread] Current Thread [Next in Thread]