- #1
Angelos K
- 48
- 0
Mod note: I revised the code below slightly, changing the loop control variable i to either j or k. The reason for this is that the browser mistakes the letter i in brackets for the BBCode italics tag, which causes some array expressions to partially disappear.
Hello,
I am trying for the first time to use LAPACK from C to diagonalize a matrix and I am stuck.
I try to modify this example http://rcabreral.blogspot.co.uk/2010/05/eigenvalues-clapack.html from zgeev to degeev. I have looked at the DGEEV input parameters, but it seems I don't understand the well enough.
Below is my code. Running it produces:
** On entry to DGEEV parameter number 9 had an illegal value
EDIT: The error occurs in the call of degeev spanning lines 48 to (including) 53.
Hello,
I am trying for the first time to use LAPACK from C to diagonalize a matrix and I am stuck.
I try to modify this example http://rcabreral.blogspot.co.uk/2010/05/eigenvalues-clapack.html from zgeev to degeev. I have looked at the DGEEV input parameters, but it seems I don't understand the well enough.
Below is my code. Running it produces:
** On entry to DGEEV parameter number 9 had an illegal value
EDIT: The error occurs in the call of degeev spanning lines 48 to (including) 53.
#include<stdio.h>
#include<math.h>
#include<complex.h>
#include <stdlib.h>
//...................
void dgeTranspose( double *Transposed, double *M ,int n)
{
int i,j;
for(i=0;i<n;i++)
for(j=0;j<n;j++) Transposed[i+n*j] = M[i*n+j];
}
//...................
// MatrixComplexEigensystem: computes the eigenvectors and eigenValues of input matrix A
// The eigenvectors are stored in columns
//..................
void MatrixComplexEigensystem( double *eigenvectorsVR, double *eigenvaluesW, double *A, int N)
{
int i;
double *AT = (double *) malloc( N*N*sizeof(double ) );
dgeTranspose( AT, A , N);
char JOBVL ='N'; // Compute Right eigenvectors
char JOBVR ='V'; // Do not compute Left eigenvectors
double VL[1];
int LDVL = 1;
int LDVR = N;
int LWORK = 4*N;
double *WORK = (double *)malloc( LWORK*sizeof(double));
double *RWORK = (double *)malloc( 2*N*sizeof(double));
int INFO;
dgeev_( &JOBVL, &JOBVR, &N, AT , &N , eigenvaluesW ,
VL, &LDVL,
eigenvectorsVR, &LDVR,
WORK,
&LWORK, RWORK, &INFO );
dgeTranspose( AT, eigenvectorsVR , N);
for(j=0;j<N*N;j++) eigenvectorsVR[j]=AT[j];
free(WORK);
free(RWORK);
free(AT);
}int main()
{
int i,j;
const int N = 3; double A[] = { 1.+I , 2. , 3 , 4. , 5.+I , 6. , 7., 8., 9. + I};
double eigenVectors[N*N];
double eigenValues[N]; MatrixComplexEigensystem( eigenVectors, eigenValues, A, N);
printf("\nEigenvectors\n");
for(k=0;k<N;k++){
for(j=0;j<N;j++) printf("%e", eigenVectors[k*N + j]);
printf("\n");
}
printf("\nEigenvalues \n");
for(j=0;i<N;j++) printf("%e", eigenValues[j] );printf("\n------------------------------------------------------------\n");
return 0;}
Any help much appreciated!#include<math.h>
#include<complex.h>
#include <stdlib.h>
//...................
void dgeTranspose( double *Transposed, double *M ,int n)
{
int i,j;
for(i=0;i<n;i++)
for(j=0;j<n;j++) Transposed[i+n*j] = M[i*n+j];
}
//...................
// MatrixComplexEigensystem: computes the eigenvectors and eigenValues of input matrix A
// The eigenvectors are stored in columns
//..................
void MatrixComplexEigensystem( double *eigenvectorsVR, double *eigenvaluesW, double *A, int N)
{
int i;
double *AT = (double *) malloc( N*N*sizeof(double ) );
dgeTranspose( AT, A , N);
char JOBVL ='N'; // Compute Right eigenvectors
char JOBVR ='V'; // Do not compute Left eigenvectors
double VL[1];
int LDVL = 1;
int LDVR = N;
int LWORK = 4*N;
double *WORK = (double *)malloc( LWORK*sizeof(double));
double *RWORK = (double *)malloc( 2*N*sizeof(double));
int INFO;
dgeev_( &JOBVL, &JOBVR, &N, AT , &N , eigenvaluesW ,
VL, &LDVL,
eigenvectorsVR, &LDVR,
WORK,
&LWORK, RWORK, &INFO );
dgeTranspose( AT, eigenvectorsVR , N);
for(j=0;j<N*N;j++) eigenvectorsVR[j]=AT[j];
free(WORK);
free(RWORK);
free(AT);
}int main()
{
int i,j;
const int N = 3; double A[] = { 1.+I , 2. , 3 , 4. , 5.+I , 6. , 7., 8., 9. + I};
double eigenVectors[N*N];
double eigenValues[N]; MatrixComplexEigensystem( eigenVectors, eigenValues, A, N);
printf("\nEigenvectors\n");
for(k=0;k<N;k++){
for(j=0;j<N;j++) printf("%e", eigenVectors[k*N + j]);
printf("\n");
}
printf("\nEigenvalues \n");
for(j=0;i<N;j++) printf("%e", eigenValues[j] );printf("\n------------------------------------------------------------\n");
return 0;}
Last edited by a moderator: