Search code examples
fortranlapackeigenvalue

Incorrect eigenvalues with cgeev

I have written the following code to get eigenvalues and determinant of the matrix in fortran with the use of cgeev:

SUBROUTINE CDETS(CS,CW,CDET,N)
IMPLICIT REAL*8 (A,B,D-H,O-Z)
IMPLICIT COMPLEX*16 (C)
DIMENSION CW(*),CS(N,*)
ALLOCATABLE :: CWK(:), WK(:), CWL(:,:),CWR(:,:)

ALLOCATE (WK(2*N),CWK(10*N),CWL(N,N),CWR(N,N))
CALL CGEEV('N','N',N,CS,N,CW,CWL,N,CWR,N,CWK,10*N,
  &         WK,INFO)

DEALLOCATE (WK,CWK,CWL,CWR)
CDET = 1.D0
DO i=1,N
  CDET = CDET*CW(i)
ENDDO
END SUBROUTINE

And this simple program to check it:

PROGRAM TESTDET
  IMPLICIT REAL*8 (A,B,D-H,O-Z)
  IMPLICIT COMPLEX*16 (C)
  DIMENSION :: CS(2,2), CW(2)
  CS(1,1)=(1.D0,1.D0)
  CS(1,2)=1.D0
  CS(2,1)=0.D0
  CS(2,2)=1.D0
  CALL CDETS(CS,CW,CDET,2)
  PRINT *, CW(1)
  PRINT *, CW(2)
END

And I get the following quite confusing results:

 (  0.0000000000000000     ,  1.0000000000000000     )
 (  1.2828908559913808E-319,  7.6130689002776223E-317)

What is the matter here?

Solution

  • You are using the wrong procedure. The argument types between cgeev and zgeev are the same but differ in kinds. For zgeev, the array RWORK is of type and kind of double precision and the arrays A, VL, VR, W, and WORK are complex*16. For cgeev the types and kinds, respectively are real and complex.

    You are implicitly typing reals as real*8, which are double precision in a non-portable fashion. Likewise your complex variables are typed as complex*16. Your variables match the argument list of zgeev, which is why it works for you and cgeev doesn't. To use cgeev change your implicit types to real and complex and you'll find cgeev now works and zgeev does not.

    From the BLAS naming conventions:

    Every routine in the BLAS library comes in four flavors, each prefixed by the letters S, D, C, and Z, respectively. Each letter indicates the format of input data:

    • S stands for single-precision (32-bit IEEE floating point numbers),
    • D stands for double-precision (64-bit IEEE floating point numbers),
    • C stands for complex numbers (represented by pairs of 32-bit IEEE floating point numbers),
    • Z stands for double complex numbers (represented by pairs of 64-bit IEEE floating point numbers)

    For your double precision complex variables, you'll need to use the Z variants of functions rather than C.