The corresponding scalapck routines to zgeev() in lapack??

[phymilton]

In my work, I need a parallel subroutine to find the eigenvalue and eigenvector of a general dense complex matrix. In serial code, I can use zgeev() to do the job, but when I parallel my current code, I cannot find a corresponding routine to zgeev().

Can you give me some advice to solve this problem?

Thank you!

[Julien Langou]

Hello,
there is no routine to compute the standard complex eigenvalue problem for a general matrix ( something like a PZGEEV ) in the current version of ScaLAPACK.
Julien Langou

[admin]

Here is a quick and dirty work around for the user who wants PZGEEV,
or indeed any missing complex version of any code.

We exploit the homomorphism between the complex number x = Rx + i*Ix where Rx is the real part of x and Ix is the imaginary part,
and the real 2-by-2 matrix

Code: Select all

 
  [ Rx  -Ix ]
  [ Ix   Rx ]


It is easy to confirm that adding, multiplying and taking reciprocals of
complex numbers corresponds to adding, multiplying and taking inverses
of their corresponding 2-by-2 matrices.

In the case of the complex non-Hermitian eigenproblem A*x = x*lambda,
we rewrite it as

Code: Select all

   [ RA  -IA ] * [ Rx  -Ix ] = [ Rx -Ix ] * [ Rlambda -Ilambda ]
   [ IA   RA ]   [ Ix   Rx ]   [ Ix  Rx ]   [ Ilambda  Rlambda ]


Since

Code: Select all

   inv(S)* [ Rlambda -Ilambda ] * S = diag( Rlambda -i* Ilambda, Rlambda + i*Ilambda )
           [ Ilambda  Rlambda ]


with S = [[1 i];[i,1]], we get

Code: Select all

   [ RA  -IA ] * [ Rx  -Ix ] * S = [ Rx -Ix ] * S *  [ conj(lambda)       0     ]
   [ IA   RA ]   [ Ix   Rx ]       [ Ix  Rx ]        [         0       lambda   ]


or

Code: Select all

   [ RA  -IA ] * [ conj(x)     i*x ] =  [ conj(x)     i*x ] * [ conj(lambda)       0     ]
   [ IA   RA ]   [ i*conj(x)     x ]    [ i*conj(x)     x ]   [         0         lambda ]


In other words, the real matrix Ahat = [[ RA, -IA];[IA, RA]] has a pair
of eigenvalues lambda and conj(lambda) corresponding to the single
eigenvalue lambda of A, with eigenvectors constructed easily from
the eigenvector x of A. So given the eigenvalues and eigenvectors
of Ahat, the eigenvalues and eigenvectors of
A can be reconstructed.

Roundoff makes all this tricky, since, for example, a real eigenvalue
of A is supposed to turn into an exact double eigenvalue of Ahat,
but roundoff may make it hard to pair them up this way.

Still, for the user this may be better than nothing.

Needless to say, this general technique can change many
complex problems into real problems, albeit by increasing
the storage and work by constant factors.

Jim

[Julien Langou]

Here is a 'quick and dirty' Lapack code to show what might be done
your ScaLapack code

Code: Select all

      PROGRAM TEST_ZGEEV_VIA_DGEEV
*    
* PZGEEV does not yet exist in the current ScaLAPACK (01/28/05), here is a     
* quick and dirty way to obatin eigencouples of a complex matrix via the real     
* routine. Demo code with Lapack.
* (Proposition of Jim Demmel)
* Looking for the eigencouples of A (complex n-by-n matrix) can be changed as
* looking for the eigncouples of
* [ Re(A) -Im(A) ]
* [ Im(A)  Re(A) ]
* see http://icl.cs.utk.edu/lapack-forum/viewtopic.php?t=18 and answers
* for more info
*
      IMPLICIT NONE 
*    
      INTEGER           N
      PARAMETER         (N=100)
      INTEGER           LWORK,LRWORK
      PARAMETER         (LWORK=2*N,LRWORK=2*N)
      INTEGER           LWORKhat
      PARAMETER         (LWORKhat=4*2*N)
*    
      INTEGER           INFO,I,J
*
      COMPLEX*16        A(N,N),VL(1,1),V(N,N),W(N),WORK(LWORK)
      COMPLEX*16        TEMP_A(N,N)
*
      DOUBLE PRECISION  Ahat(2*N,2*N),WRhat(2*N),WIhat(2*N),VLhat(1,1)
      DOUBLE PRECISION  Vhat(2*N,2*N), WORKhat(LWORKhat)
      DOUBLE PRECISION  TEMP_Ahat(2*N,2*N)
      COMPLEX*16        VV1(N,N),VV2(N,N),AVV1(N,N),AVV2(N,N)
      COMPLEX*16        VV(N,N),WW(N)
*
      DOUBLE PRECISION  RWORK(LRWORK),ERROR,ERROR_
      DOUBLE PRECISION  ERROR1,ERROR2,NORMVV
      COMPLEX*16        ZMONE,ZZERO,LAMBDA1,LAMBDA2,ZONE
*   
      EXTERNAL          DGEEV,ZGEEV
*    
      COMPLEX*16        ZDOTC
      DOUBLE PRECISION  DZNRM2
      EXTERNAL          DZNRM2,ZDOTC
*    
      INTRINSIC         RAND,ABS,DCMPLX,DBLE,DIMAG
*    
      WRITE(*,*) 'ZGEEV via DGEEV, N=',N
*
      ZMONE = DCMPLX(-1.0d0,0.0d0)
      ZONE =  DCMPLX( 1.0d0,0.0d0)
      ZZERO = DCMPLX( 0.0d0,0.0d0)
*    
      DO J=1,N
         DO I=1,N
            A(I,J)  = DCMPLX(RAND(),RAND())
         END DO   
      END DO   
*
      DO J=1,N
         DO I=1,N
            Ahat(I,J)     = DBLE(A(I,J))
         END DO   
      END DO   
      DO J=1,N
         DO I=1,N
            Ahat(N+I,N+J) = (Ahat(I,J))
         END DO   
      END DO   
      DO J=1,N
         DO I=1,N
            Ahat(N+I,J)   = DIMAG(A(I,J))
         END DO   
      END DO   
      DO J=1,N
         DO I=1,N
            Ahat(I,N+J)   = -Ahat(N+I,J)
         END DO   
      END DO   
*    
      CALL ZCOPY(2*N*2*N,Ahat,1,TEMP_Ahat,1)
      CALL DGEEV('N','V',2*N,TEMP_Ahat,2*N,
     &           WRhat,WIhat,VLhat,1,Vhat,2*N,
     &           WORKhat,LWORKhat,INFO)
*
      DO J=1,N
         DO I=1,N
            VV1(I,J) = DCMPLX(Vhat(I,2*J-1),Vhat(I,2*J))
         END DO
      END DO
      DO J=1,N
         DO I=1,N
            VV2(I,J) = DCMPLX(Vhat(I,2*J-1),-Vhat(I,2*J))
         END DO
      END DO
*
      CALL ZGEMM('N','N',N,N,N,ZONE,A,N,VV1,N,ZZERO,AVV1,N)
      CALL ZGEMM('N','N',N,N,N,ZONE,A,N,VV2,N,ZZERO,AVV2,N)
*
      ERROR = 0.d0
      DO J=1,N
         NORMVV  = DZNRM2(N,VV2(1,J),1)
         LAMBDA1 = ZZERO
         DO I=1,N
            LAMBDA1 = LAMBDA1 + CONJG(VV1(I,J))*AVV1(I,J)
         END DO
*         LAMBDA1 = ZDOTC(N,VV1(1,J),1,AVV1(1,J),1)
         LAMBDA1 = LAMBDA1/NORMVV/NORMVV
         CALL ZAXPY(N,-LAMBDA1,VV1(1,J),1,AVV1(1,J),1)
         ERROR1 = DZNRM2(N,AVV1(1,J),1)
         ERROR1 = ERROR1/CDABS(LAMBDA1)/NORMVV
*
         LAMBDA2 = ZZERO
         DO I=1,N
            LAMBDA2 = LAMBDA2 + CONJG(VV2(I,J))*AVV2(I,J)
         END DO
*         LAMBDA2 = ZDOTC(N,VV2(1,J),1,AVV2(1,J),1)
         LAMBDA2 = LAMBDA2/NORMVV/NORMVV
         CALL ZAXPY(N,-LAMBDA2,VV2(1,J),1,AVV2(1,J),1)
         ERROR2 = DZNRM2(N,AVV2(1,J),1)
         ERROR2 = ERROR2/CDABS(LAMBDA2)/NORMVV
         IF (ERROR1>ERROR2) THEN
            CALL ZCOPY(N,VV2(1,J),1,VV(1,J),1)
            WW(J) = LAMBDA2
            IF (ERROR2>ERROR) THEN
              ERROR = ERROR2
            END IF
         ELSE
            CALL ZCOPY(N,VV1(1,J),1,VV(1,J),1)
            WW(J) = LAMBDA1
            IF (ERROR1>ERROR) THEN
              ERROR = ERROR1
            END IF
         END IF
      END DO   
      WRITE(*,*) 'max relative bwd error',ERROR
*
      WRITE(*,*) 'standard ZGEEV, N=',N
      CALL ZCOPY(N*N,A,1,TEMP_A,1)
      CALL ZGEEV('N','V',N,TEMP_A,N,W,VL,1,V,N,WORK,LWORK,RWORK,INFO)
*
      CALL ZGEMM('N','N',N,N,N,ZMONE,A,N,V,N,ZZERO,TEMP_A,N)
*
      ERROR = 0.d0
      DO J=1,N
         CALL ZAXPY(N,W(J),V(1,J),1,TEMP_A(1,J),1)
         ERROR_ = DZNRM2(N,TEMP_A(1,J),1)/CDABS(W(J))/DZNRM2(N,V(1,J),1)
         IF (ERROR_>ERROR) THEN
            ERROR = ERROR_   
         END IF
      END DO   
      WRITE(*,*) 'max relative bwd error',ERROR
*
      END 


[phymilton]

Thank you very much for helping on this topic!

I also search the Internet for result: here is one and the only one I found that provide some info about the non-symmetric complex eigen solver:

Can you help me to justify that which method is better to do the job: (power point page 39)

http://www.hpc.unm.edu/news_events_workshops/workshops/parallel_numerical_libraries/slides/libs.ppt

Thanks again.

Ming

[phymilton]

I checked all the routines, there is no pdgeev either. Is this mean, I have to go from all the building blocks of scalapck to make a pdgeev function routine then convert it to a pzgeev function? I think I cannot go through all this according to my knowledge of programming and linear algebra.

Does the method provide in the link work for large arraies? As I test the sample programs download form http://www.ncsa.uiuc.edu's webct course, it can only run for small array. When the array goes larger, it crashes.

Thank you again!

[Julien Langou]

Is this mean, I have to go from all the building blocks of scalapack to make a pdgeev function routine then convert it to a pzgeev function?

Yes that's the best we can propose with the current version of ScaLAPACK.
Neither the routine pzgeev does exist, nor the routine pdgeev. But the real eigenvalue problem can be tackled out with ScaLAPACK in two steps:
- reduction to upper Hessenberg form by PDGEHRD, followed by,
- reduction of the Hessenberg matrix to Schur form by PDLAHQR.
Then the complex eigenvalue problem can be tacked out from real eigenvalue problem as explained above.
Julien Langou