Using expm in C++

classic Classic list List threaded Threaded
14 messages Options
Reply | Threaded
Open this post in threaded view
|

Using expm in C++

loisp
Hello,

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;
                }
        }

        hamiltonian = hamiltonian.expm();


When I try to compile my program using mkoctfile I get the following error

        dmatrix.cpp: In function ‘int main()’:
        dmatrix.cpp:54: error: ‘class Matrix’ has no member named ‘expm’

I was wondering if the classes Matrix and ComplexMatrix provide that function and if they don't, do you have any work around to suggest? I am using MacOSX and have Octave installed though fink.

Thanks in advance.

Lois
Reply | Threaded
Open this post in threaded view
|

Re: Using expm in C++

loisp
This post was updated on .
This is the link to the octave package database in Fink.

I have octave360, octave360-dev, octave360-shlibs, linear-algebra-oct360, gsl-oct360 installed. My code includes

#include <octave/oct.h>
#include <octave/dMatrix.h> 
#include <octave/CMatrix.h>

Reply | Threaded
Open this post in threaded view
|

Re: Using expm in C++

Alexander Hansen-2
-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

On 2/6/12 6:06 PM, loisp wrote:
> http://pdb.finkproject.org/pdb/browse.php?summary=octave This is
> the link to the octave package database in Fink.
>
> I have octave360, octave360-dev, octave360-shlibs,
> linear-algebra-oct360, gsl-oct360 installed.
>
> --

The general Octave populace my not know what these mean. :-)

The important part is that "octave360-dev" contains all of the
headers, and I made no modifications to them.

There's an "expm" function, but it's implemented as an m-file rather
than in the C++ API, as far as I can tell--definitely not in Matrix.h:

$ dpkg -L octave360-atlas-dev | xargs grep expm
/sw/include/octave-3.6.0/octave/config.h:/* Define to 1 if you have
the `expm1' function. */
/sw/include/octave-3.6.0/octave/config.h:/* Define to 1 if you have
the `expm1f' function. */
/sw/include/octave-3.6.0/octave/lo-specfun.h:extern OCTAVE_API double
expm1 (double x);
/sw/include/octave-3.6.0/octave/lo-specfun.h:extern OCTAVE_API Complex
expm1 (const Complex& x);
/sw/include/octave-3.6.0/octave/lo-specfun.h:extern OCTAVE_API float
expm1f (float x);
/sw/include/octave-3.6.0/octave/lo-specfun.h:extern OCTAVE_API
FloatComplex expm1 (const FloatComplex& x);
/sw/include/octave-3.6.0/octave/ov-base.h:      umap_expm1,
/sw/include/octave-3.6.0/octave/ov.h:  MAPPER_FORWARD (expm1)
/sw/share/doc/octave360-atlas-dev/ChangeLog:
/linear-algebra/duplication_matrix.m scripts/linear-algebra/expm.m
/sw/share/doc/octave360-atlas-dev/NEWS:      expm1
nargoutchk                  strchr
- --
Alexander Hansen, Ph.D.
Fink User Liaison
http://finkakh.wordpress.com/
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.11 (Darwin)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/

iEYEARECAAYFAk8wcUIACgkQB8UpO3rKjQ+fDACfVTnXVbkDyS6rE0XNrJOPBrk6
FUcAniLAe+8SlcDiMdwcHawIiFPwnPWj
=ZfiH
-----END PGP SIGNATURE-----
_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: Using expm in C++

loisp
Yes, there is just an m-file. Do you have any suggestions how I can use the expm function in my C++ code?

On Tue, Feb 7, 2012 at 6:04 AM, Alexander Hansen-2 [via Octave] <[hidden email]> wrote:
-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

On 2/6/12 6:06 PM, loisp wrote:
> http://pdb.finkproject.org/pdb/browse.php?summary=octave This is
> the link to the octave package database in Fink.
>
> I have octave360, octave360-dev, octave360-shlibs,
> linear-algebra-oct360, gsl-oct360 installed.
>
> --

The general Octave populace my not know what these mean. :-)

The important part is that "octave360-dev" contains all of the
headers, and I made no modifications to them.

There's an "expm" function, but it's implemented as an m-file rather
than in the C++ API, as far as I can tell--definitely not in Matrix.h:

$ dpkg -L octave360-atlas-dev | xargs grep expm
/sw/include/octave-3.6.0/octave/config.h:/* Define to 1 if you have
the `expm1' function. */
/sw/include/octave-3.6.0/octave/config.h:/* Define to 1 if you have
the `expm1f' function. */
/sw/include/octave-3.6.0/octave/lo-specfun.h:extern OCTAVE_API double
expm1 (double x);
/sw/include/octave-3.6.0/octave/lo-specfun.h:extern OCTAVE_API Complex
expm1 (const Complex& x);
/sw/include/octave-3.6.0/octave/lo-specfun.h:extern OCTAVE_API float
expm1f (float x);
/sw/include/octave-3.6.0/octave/lo-specfun.h:extern OCTAVE_API
FloatComplex expm1 (const FloatComplex& x);
/sw/include/octave-3.6.0/octave/ov-base.h:      umap_expm1,
/sw/include/octave-3.6.0/octave/ov.h:  MAPPER_FORWARD (expm1)
/sw/share/doc/octave360-atlas-dev/ChangeLog:
/linear-algebra/duplication_matrix.m scripts/linear-algebra/expm.m
/sw/share/doc/octave360-atlas-dev/NEWS:      expm1
nargoutchk                  strchr
- --
Alexander Hansen, Ph.D.
Fink User Liaison
http://finkakh.wordpress.com/
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.11 (Darwin)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/

iEYEARECAAYFAk8wcUIACgkQB8UpO3rKjQ+fDACfVTnXVbkDyS6rE0XNrJOPBrk6
FUcAniLAe+8SlcDiMdwcHawIiFPwnPWj
=ZfiH
-----END PGP SIGNATURE-----
_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave



If you reply to this email, your message will be added to the discussion below:
http://octave.1599824.n4.nabble.com/Using-expm-in-C-tp4363138p4363379.html
To unsubscribe from Using expm in C++, click here.
NAML

Reply | Threaded
Open this post in threaded view
|

Re: Using expm in C++

loisp
In reply to this post by Alexander Hansen-2
Btw thanks for your reply :)

On Tue, Feb 7, 2012 at 6:03 AM, Alexander Hansen <[hidden email]> wrote:
-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

On 2/6/12 6:06 PM, loisp wrote:
> http://pdb.finkproject.org/pdb/browse.php?summary=octave This is
> the link to the octave package database in Fink.
>
> I have octave360, octave360-dev, octave360-shlibs,
> linear-algebra-oct360, gsl-oct360 installed.
>
> --

The general Octave populace my not know what these mean. :-)

The important part is that "octave360-dev" contains all of the
headers, and I made no modifications to them.

There's an "expm" function, but it's implemented as an m-file rather
than in the C++ API, as far as I can tell--definitely not in Matrix.h:

$ dpkg -L octave360-atlas-dev | xargs grep expm
/sw/include/octave-3.6.0/octave/config.h:/* Define to 1 if you have
the `expm1' function. */
/sw/include/octave-3.6.0/octave/config.h:/* Define to 1 if you have
the `expm1f' function. */
/sw/include/octave-3.6.0/octave/lo-specfun.h:extern OCTAVE_API double
expm1 (double x);
/sw/include/octave-3.6.0/octave/lo-specfun.h:extern OCTAVE_API Complex
expm1 (const Complex& x);
/sw/include/octave-3.6.0/octave/lo-specfun.h:extern OCTAVE_API float
expm1f (float x);
/sw/include/octave-3.6.0/octave/lo-specfun.h:extern OCTAVE_API
FloatComplex expm1 (const FloatComplex& x);
/sw/include/octave-3.6.0/octave/ov-base.h:      umap_expm1,
/sw/include/octave-3.6.0/octave/ov.h:  MAPPER_FORWARD (expm1)
/sw/share/doc/octave360-atlas-dev/ChangeLog:
/linear-algebra/duplication_matrix.m scripts/linear-algebra/expm.m
/sw/share/doc/octave360-atlas-dev/NEWS:      expm1
nargoutchk                  strchr
- --
Alexander Hansen, Ph.D.
Fink User Liaison
http://finkakh.wordpress.com/
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.11 (Darwin)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/

iEYEARECAAYFAk8wcUIACgkQB8UpO3rKjQ+fDACfVTnXVbkDyS6rE0XNrJOPBrk6
FUcAniLAe+8SlcDiMdwcHawIiFPwnPWj
=ZfiH
-----END PGP SIGNATURE-----


_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: Using expm in C++

bpabbott
Administrator
In reply to this post by loisp

On Feb 6, 2012, at 5:58 PM, loisp wrote:

> Hello,
>
> 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;
> }
> }
>
> hamiltonian = hamiltonian.expm();
>
>
> When I try to compile my program using mkoctfile I get the following error
>
> dmatrix.cpp: In function ‘int main()’:
> dmatrix.cpp:54: error: ‘class Matrix’ has no member named ‘expm’
>
> I was wondering if the classes Matrix and ComplexMatrix provide that
> function and if they don't, do you have any work around to suggest? I am
> using MacOSX and have Octave installed though fink.
>
> Thanks in advance.
>
> Lois

I may be wrong, but I don't think there is an expm is an on the c++ side. However, there is an m-file version. You can call it via feval (although I've never tried).

http://www.gnu.org/software/octave/doc/interpreter/Calling-Octave-Functions-from-Oct_002dFiles.html

Ben


_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: Using expm in C++

loisp
The thing is I don't want to use an oct-file. I want my code to be in C++ and to use Octave just for the expm() function. Is there any way I can achieve that?

On Tue, Feb 7, 2012 at 6:48 AM, bpabbott [via Octave] <[hidden email]> wrote:

On Feb 6, 2012, at 5:58 PM, loisp wrote:

> Hello,
>
> 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;
> }
> }
>
> hamiltonian = hamiltonian.expm();
>
>
> When I try to compile my program using mkoctfile I get the following error
>
> dmatrix.cpp: In function ‘int main()’:
> dmatrix.cpp:54: error: ‘class Matrix’ has no member named ‘expm’
>
> I was wondering if the classes Matrix and ComplexMatrix provide that
> function and if they don't, do you have any work around to suggest? I am
> using MacOSX and have Octave installed though fink.
>
> Thanks in advance.
>
> Lois
I may be wrong, but I don't think there is an expm is an on the c++ side. However, there is an m-file version. You can call it via feval (although I've never tried).

http://www.gnu.org/software/octave/doc/interpreter/Calling-Octave-Functions-from-Oct_002dFiles.html

Ben


_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave



If you reply to this email, your message will be added to the discussion below:
http://octave.1599824.n4.nabble.com/Using-expm-in-C-tp4363138p4363437.html
To unsubscribe from Using expm in C++, click here.
NAML

Reply | Threaded
Open this post in threaded view
|

Re: Using expm in C++

bpabbott
Administrator

On Feb 6, 2012, at 8:24 PM, loisp wrote:

> On Tue, Feb 7, 2012 at 6:48 AM, bpabbott [via Octave] <[hidden email]> wrote:
>
> On Feb 6, 2012, at 5:58 PM, loisp wrote:
>
> > Hello,
> >
> > 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;
> > }
> > }
> >
> > hamiltonian = hamiltonian.expm();
> >
> >
> > When I try to compile my program using mkoctfile I get the following error
> >
> > dmatrix.cpp: In function ‘int main()’:
> > dmatrix.cpp:54: error: ‘class Matrix’ has no member named ‘expm’
> >
> > I was wondering if the classes Matrix and ComplexMatrix provide that
> > function and if they don't, do you have any work around to suggest? I am
> > using MacOSX and have Octave installed though fink.
> >
> > Thanks in advance.
> >
> > Lois
> I may be wrong, but I don't think there is an expm is an on the c++ side. However, there is an m-file version. You can call it via feval (although I've never tried).
>
> http://www.gnu.org/software/octave/doc/interpreter/Calling-Octave-Functions-from-Oct_002dFiles.html
>
> Ben
>
> The thing is I don't want to use an oct-file. I want my code to be in C++ and to use Octave just for the expm() function. Is there any way I can achieve that?

I'm unfamiliar with what you're asking, but I'll try to help.

If you only need to call expm and don't require anything else, I think you're better off recoding the expm function in c++.

        http://hg.savannah.gnu.org/hgweb/octave/file/3e4350f09a55/scripts/linear-algebra/expm.m

If you really want to call Octave you can start with ...

        http://octave.sourceforge.net/doxygen/html/

There's a rather dated description for call Octave from c++ below.

        http://www.mathias-michel.de/download/howto-octave-c++.ps

Ben

_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: Using expm in C++

loisp
This post was updated on .
I need to call just expm() and I don't know how to write a C++ API for
expm.m. I have also tried Mathias Michel's guide without any success.

Mine is a standalone program. I have also found an old thread where the same issue was raised.

On Tue, Feb 7, 2012 at 7:18 AM, bpabbott [via Octave] <
ml-node+s1599824n4363490h41@n4.nabble.com> wrote:

>
> On Feb 6, 2012, at 8:24 PM, loisp wrote:
>
> > On Tue, Feb 7, 2012 at 6:48 AM, bpabbott [via Octave] <[hidden email]>
> wrote:
> >
> > On Feb 6, 2012, at 5:58 PM, loisp wrote:
> >
> > > Hello,
> > >
> > > 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;
> > > }
> > > }
> > >
> > > hamiltonian = hamiltonian.expm();
> > >
> > >
> > > When I try to compile my program using mkoctfile I get the following
> error
> > >
> > > dmatrix.cpp: In function ‘int main()’:
> > > dmatrix.cpp:54: error: ‘class Matrix’ has no member named ‘expm’
> > >
> > > I was wondering if the classes Matrix and ComplexMatrix provide that
> > > function and if they don't, do you have any work around to suggest? I
> am
> > > using MacOSX and have Octave installed though fink.
> > >
> > > Thanks in advance.
> > >
> > > Lois
> > I may be wrong, but I don't think there is an expm is an on the c++
> side. However, there is an m-file version. You can call it via feval
> (although I've never tried).
> >
> >
> http://www.gnu.org/software/octave/doc/interpreter/Calling-Octave-Functions-from-Oct_002dFiles.html
> >
> > Ben
> >
> > The thing is I don't want to use an oct-file. I want my code to be in
> C++ and to use Octave just for the expm() function. Is there any way I can
> achieve that?
>
> I'm unfamiliar with what you're asking, but I'll try to help.
>
> If you only need to call expm and don't require anything else, I think
> you're better off recoding the expm function in c++.
>
>
> http://hg.savannah.gnu.org/hgweb/octave/file/3e4350f09a55/scripts/linear-algebra/expm.m
>
> If you really want to call Octave you can start with ...
>
>         http://octave.sourceforge.net/doxygen/html/
>
> There's a rather dated description for call Octave from c++ below.
>
>         http://www.mathias-michel.de/download/howto-octave-c++.ps
>
> Ben
>
> _______________________________________________
> Help-octave mailing list
> [hidden email] <http://user/SendEmail.jtp?type=node&node=4363490&i=0>
> https://mailman.cae.wisc.edu/listinfo/help-octave
>
>
> ------------------------------
>  If you reply to this email, your message will be added to the discussion
> below:
> http://octave.1599824.n4.nabble.com/Using-expm-in-C-tp4363138p4363490.html
>  To unsubscribe from Using expm in C++, click here<http://octave.1599824.n4.nabble.com/template/NamlServlet.jtp?macro=unsubscribe_by_code&node=4363138&code=bG9pcy5zaGVraW5haEBnbWFpbC5jb218NDM2MzEzOHw0MDY5MzcyODM=>
> .
> NAML<http://octave.1599824.n4.nabble.com/template/NamlServlet.jtp?macro=macro_viewer&id=instant_html%21nabble%3Aemail.naml&base=nabble.naml.namespaces.BasicNamespace-nabble.view.web.template.NabbleNamespace-nabble.view.web.template.NodeNamespace&breadcrumbs=notify_subscribers%21nabble%3Aemail.naml-instant_emails%21nabble%3Aemail.naml-send_instant_email%21nabble%3Aemail.naml>
>
Reply | Threaded
Open this post in threaded view
|

Re: Using expm in C++

sustik
In reply to this post by loisp
Hi,

Are you sure you need full matrix exponentiation? That is for arbitrary matrices?  If your matrices were symmetric then I would do an eigendecomposition VDV' and compute Vexp(D)V'.  When I need to do exponentiation many times in the algorithm I try to maintain the eigendecomposition throughout.

I have done this in my research in machine learning applications.  If you tell me more about your overall problem I can try to help further. You can ommit the list.

Good luck!

-Matyas 


Sent from an Android phone running a Linux kernel.

_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: Using expm in C++

c.-2
In reply to this post by loisp
On 7 Feb 2012, at 03:07, loisp wrote:

> I need to call just expm() and I don't know how to write a C++ API for expm.m. I have also tried Mathias Michel's guide without any success.

If I understand correctly what you want to do is write a stand-alone C++ program that can call functions defined as m-files.
What you want to do is usually referred to as "embedding" the Octave interpreter in your C++ application, if you google that phrase
you'll find more info and examples, but to get started you can just have a look at this example in the Octave source tree:

http://hg.savannah.gnu.org/hgweb/octave/file/3e4350f09a55/examples/embedded.cc

the essential parts of the example are:

1) the call to "octave_main" on line 13:

  octave_main (2, argv.c_str_vec(), 1);

 which initializes the Octave intepreter

2) lines 29 and 30 that set up an "octave_value_list" object with the inputs for the m-file function,
   call the function via "feval" and collect the output in another "octave_value_list" object.

Beware that using this feature the performance of the code will be limited by the interpreter speed,
so if you use this a lot you might be better off writing your whole code in the Octave interpreter language
rather than in C++ as using the C++ API will not give you any speed-up.

HTH,
c.


_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: Using expm in C++

loisp
In reply to this post by sustik
Thanks! That did it

Reply | Threaded
Open this post in threaded view
|

Re: Using expm in C++

sustik
In reply to this post by loisp
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"
#endif

// 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
    startgdb();
#endif

    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...
    free(iwork);
    free(work);
    free(isuppz);
    free(e);
    free(d);
   
    // 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;
    }
    free(w);

    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);

    free(z1);
    free(z);
    return;
}

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]));
    return;
}

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)

CXX=g++

CXXFLAGS = -Wall -fpic -pthread -shared  -fno-omit-frame-pointer -ansi -D_GNU_SOURCE -D_FILE_OFFSET_BITS=64
CXXOPTFLAGS = -O3 $(CXXFLAGS)
CXXDBGFLAGS = -g -DGDEBUG $(CXXFLAGS)

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

.SUFFIXES:

# 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 $@

.SECONDARY :

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,
-Matyas

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]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: Using expm in C++

loisp
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++) 

   { 

   if(i==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). 

Cheers,
Lois

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"
#endif

// 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
    startgdb();
#endif

    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...
    free(iwork);
    free(work);
    free(isuppz);
    free(e);
    free(d);
   
    // 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;
    }
    free(w);

    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);

    free(z1);
    free(z);
    return;
}

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]));
    return;
}

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)

CXX=g++

CXXFLAGS = -Wall -fpic -pthread -shared  -fno-omit-frame-pointer -ansi -D_GNU_SOURCE -D_FILE_OFFSET_BITS=64
CXXOPTFLAGS = -O3 $(CXXFLAGS)
CXXDBGFLAGS = -g -DGDEBUG $(CXXFLAGS)

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

.SUFFIXES:

# 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 $@

.SECONDARY :

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,
-Matyas

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]
https://mailman.cae.wisc.edu/listinfo/help-octave



If you reply to this email, your message will be added to the discussion below:
http://octave.1599824.n4.nabble.com/Using-expm-in-C-tp4363138p4380313.html
To unsubscribe from Using expm in C++, click here.
NAML