how to get eigenvalues of symetric matrix

[bruce]

Hi,

I want to get a eigenvalues of symetric matrix, but I have no idea which routine i should use.

I want to call the function in C code to get the result.

Any help will be highly appreciated.

[Julie]

Have a look at the LAPACK routines:
dsyev (for double, ssyev if you are using float).
They are also more advanced routines:
dysevr, dsyevd, dsyevx. You can have a look at them as well. You should have no problem to call the LAPACK library from C, just keep in mind that LAPACK is a Fortran library.
Julie

[icon]

Have a look at the LAPACK routines:
dsyev (for double, ssyev if you are using float).
They are also more advanced routines:
dysevr, dsyevd, dsyevx. You can have a look at them as well. You should have no problem to call the LAPACK library from C, just keep in mind that LAPACK is a Fortran library.
Julie

How exactly do you store A? I don't understand this:

"On entry, the symmetric matrix A. If UPLO = 'U',
the leading N-by-N upper triangular part of A con-
tains the upper triangular part of the matrix A.
If UPLO = 'L', the leading N-by-N lower triangular
part of A contains the lower triangular part of
the matrix A. On exit, if JOBZ = 'V', then if
INFO = 0, A contains the orthonormal eigenvectors
of the matrix A. If JOBZ = 'N', then on exit the
lower triangle (if UPLO='L') or the upper triangle
(if UPLO='U') of A, including the diagonal, is
destroyed.

"

[Julie]

You need to input A in column major format. This means you describe A column by column. A needs to be of size at least N^2. Since A is symmetric, you only need to fill either the upper part or the lower part. (Well if you fill the whole matrix it is not a big deal). If you have filled the lower part of the matrix you set UPLO to 'L'. If you have filled the upper part then you set UPLO to 'L'.

A quick and dirty C code can be:

Code: Select all

#include <stdio.h>
#include <stdlib.h>

extern void dsyev_( char *jobz, char *uplo, int *n, double *a, int *lda,
        double *w, double *work, int *lwork, int *info );


int main(){

        double *a=NULL, *work=NULL, *w=NULL;
        int n,lwork,info;

        n = 3;
        a = (double *)malloc(n*n*sizeof(double));
        w = (double *)malloc(n*sizeof(double));

        a[0] = 2;
        a[1] = 1; a[4] = 2;
        a[2] = 1; a[5] = 1; a[8] = 2;

        work = (double *)malloc(1*sizeof(double));
        lwork=-1;
        dsyev_( "N", "L", &n, a, &n, w, work, &lwork, &info );
        lwork= work[0];
        free(work);
        work = (double *)malloc(lwork*sizeof(double));

        dsyev_( "N", "L", &n, a, &n, w, work, &lwork, &info );

        fprintf(stdout,"e1=%f\n",w[0]);
        fprintf(stdout,"e2=%f\n",w[1]);
        fprintf(stdout,"e3=%f\n",w[2]);

        free(work);
        free(w);
        free(a);

        return 0;

}



it compute the eigenvalue of the matrix

Code: Select all

A = [ 2 1 1
      1 2 1
      1 1 2
]


Answer is 1,1,4.
Julie.

[bruce]

I got it. I will try.

[icon]

You need to input A in column major format. This means you describe A column by column. A needs to be of size at least N^2. Since A is symmetric, you only need to fill either the upper part or the lower part. (Well if you fill the whole matrix it is not a big deal). If you have filled the lower part of the matrix you set UPLO to 'L'. If you have filled the upper part then you set UPLO to 'L'.
Julie.

Cool. I see, so A can not be the original symmetric matrix?

Or can it be since you said "(Well if you fill the whole matrix it is not a big deal)"? If I use the original then 'U' or 'L' does not apply.....?????

Thank you!

[Julie]

You can use the original matrix if you want then you pick UPLO = 'U' or 'V'. The one you prefer. If you pick U, DSYEV will only consider the upper part of the matrix and undertsand that the lower part is the same. If you pick L, DSYEV will only consider the lower part of the matrix and understand the upper part is the same.
Julie.

[benndamann33]

I am having trouble getting lapack packages to work with C++ at all. Where should I have lapack downloaded to in the first place? The message that I receive is the foolowing:
internet.cpp:41:3: warning: no newline at end of file
internet.cpp: In function ‘int main()’:
internet.cpp:25: warning: converting to ‘int’ from ‘double’
/tmp/ccUGZQ3o.o: In function `main':internet.cpp:(.text+0xfc): undefined reference to `dsyev_(char*, char*, int*, double*, int*, double*, double*, int*, int*)'
:internet.cpp:(.text+0x17e): undefined reference to `dsyev_(char*, char*, int*, double*, int*, double*, double*, int*, int*)'
collect2: ld returned 1 exit status

where I have named the example file given in this forum internet.cpp

[Julie]

Hello:

oups sorry, this code was meant for C not C++.
The main error:

Code: Select all

/tmp/ccUGZQ3o.o: In function `main':internet.cpp:(.text+0xfc): undefined reference to `dsyev_(char*, char*, int*, double*, int*, double*, double*, int*, int*)'
:internet.cpp:(.text+0x17e): undefined reference to `dsyev_(char*, char*, int*, double*, int*, double*, double*, int*, int*)'
collect2: ld returned 1 exit status


comes from the fact that in C++ you need to add "C" when declaring the C or Fortran functions as external. So lines 4 and 5 should be:

Code: Select all

extern void dsyev_( char *jobz, char *uplo, int *n, double *a, int *lda,
       double *w, double *work, int *lwork, int *info );


You can remove the warning

Code: Select all

internet.cpp:25: warning: converting to ‘int’ from ‘double’


by adding a cast at line 25, this gives:

Code: Select all

        lwork= (int) work[0];
[

This is cleaner this way, yes.

I do not really understand the warning:

Code: Select all

internet.cpp:41:3: warning: no newline at end of file


? may be remove the last blank line ?

Finally, you need to link with your LAPACK and BLAS library, and if your linker is the C++ linker, then you will also need to link with the Fortran runtime libraries (for example libg2c if you are using g77 to compile LAPACK) and the math library ( -lm ).

Gives something like:

Code: Select all

g++  example_dsyev.cpp -o  ${LAPACK_DIR}/liblapack.a ${BLASDIR}/libblas.a  -lg2c -lm


Julie

[benndamann33]

thanks a lot. That warning about the line at the end of the file...should have removed that before pasting, my compiler always gives me that it's no big deal. Thanks again
Ben

[benndamann33]

oh actually..i'm a little confused


comes from the fact that in C++ you need to add "C" when declaring the C or Fortran functions as external. So lines 4 and 5 should be:
Code:

extern void dsyev_( char *jobz, char *uplo, int *n, double *a, int *lda,
double *w, double *work, int *lwork, int *info );


I don't see where the C was added in that code?

[benndamann33]

I changed it to extern "C" instead of just extern, becuase I think that's what you meant and i was being dumb. However, i'm now getting the following erros which I do not understand

g++ internet.cpp -o a.out /usr/lib/liblapack.a /usr/lib/libblas.a -lg2c -lm
/usr/lib/libblas.a(xerbla.o): In function `xerbla_': multiple definition of `xerbla_'
/usr/lib/liblapack.a(xerbla.o): first defined here
/usr/bin/ld: Warning: size of symbol `xerbla_' changed from 73 in /usr/lib/liblapack.a(xerbla.o) to 55 in /usr/lib/libblas.a(xerbla.o)
collect2: ld returned 1 exit status

[Julie]

Good interpretation. We are almost there. So now, your error message is nastier. Since you are using g++, I aussme you have also used g77. Can you try to link with the 77 linker rather than the ++:

Code: Select all

g77 internet.cpp -o a.out /usr/lib/liblapack.a /usr/lib/libblas.a


???

[benndamann33]

following error...I am not sure how to fix because I simply apt-get installed g77...will I have ot download it and edit the make file to get this to work?


g77: installation problem, cannot exec `cc1plus': No such file or directory

[Julie]

Could you tell me which lapack libray and which blas library you are using?
Is it the lapack in Fortran77 that you download from netlib or something else?
Is the blas the reference blas that is included in the LAPACK distribution or is it a vendor or optimized blas for your machine? Is it in Fortran, C or C++?

Julie