[Top][All Lists]

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

Re: Using expm in C++

From: Matyas Sustik
Subject: Re: Using expm in C++
Date: Sat, 11 Feb 2012 18:31:57 -0600

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.

reply via email to

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