On entry to DGEEV parameter number 9 had an illegal value

I am trying for the first time to use LAPACK from C to diagonalize a matrix and I am stuck.

I have been trying to modify this example http://rcabreral.blogspot.co.uk/2010/05/eigenvalues-clapack.html from zgeev to dgeev. I have looked at the DGEEV input parameters, http://www.netlib.org/lapack/explore-html/d9/d28/dgeev_8f.html but it seems I don't understand the well enough.

Hence, the code below produces:

**** On entry to DGEEV parameter number 9 had an illegal value**

EDIT: The error occurs in the call of dgeev spanning lines 48 to (including) 53.

EDIT: Note that the arguments differ from the specifications here http://www.netlib.org/lapack/explore-html/d9/d28/dgeev_8f.html in that they have been translated to pointers. That is necessary when using these Fortran routines in C, as explained here: http://www.physics.orst.edu/~rubin/nacphy/lapack/cprogp.html

#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(i=0;i<N*N;i++) eigenvectorsVR[i]=AT[i];

  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(i=0;i<N;i++){
    for(j=0;j<N;j++) printf("%e", eigenVectors[i*N + j]);
    printf("\n");
  }
  printf("\nEigenvalues \n");
  for(i=0;i<N;i++) printf("%e",  eigenValues[i] );
  printf("\n------------------------------------------------------------\n"); 

  return 0;
}

Solution

You can not port directly from zgeev to dgeev. The zgeev gets a complex matrix and computes complex eigenvalues. While dgeev gets a real matrix and computes complex eigenvalues. In order to be consistent LAPACK uses WR and WI which is used for the real and imaginary part of each eigenvalue.

So note that dgeev definition is

void dgeev_(char* JOBVL, char* JOBVR, int* N, double* A, int* LDA, double* WR, double* WI, double* VL, int* LDVL, double* VR, int* LDVR, double* WORK, int* LWORK, int* INFO);

My suggestion for your example is to remove:

#include <complex.h>

remove I's from matrix of doubles:

double A[] = { 1. , 2. ,  3 , 4. , 5. , 6. , 7., 8., 9.};

then double the size of eigenvalues vector:

double eigenValues[2*N];

and call dgeev using WR and WI:

double *eigenvaluesWR = eigenvaluesW;
double *eigenvaluesWI = eigenvaluesW+N;
dgeev_(&JOBVL, &JOBVR, &N, AT, &N, 
    eigenvaluesWR, eigenvaluesWI, 
    VL, &LDVL, 
    eigenvectorsVR, &LDVR,
    WORK, &LWORK, &INFO);