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

Open discussion regarding features, bugs, issues, vendors, etc.

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

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!
phymilton

Posts: 19
Joined: Mon Jan 24, 2005 11:41 pm
Location: Ames, IA

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
Julien Langou

Posts: 834
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

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

Posts: 616
Joined: Wed Dec 08, 2004 7:07 pm

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

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  `
Julien Langou

Posts: 834
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

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

Posts: 19
Joined: Mon Jan 24, 2005 11:41 pm
Location: Ames, IA

### One more problem

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!
phymilton

Posts: 19
Joined: Mon Jan 24, 2005 11:41 pm
Location: Ames, IA

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
Julien Langou

Posts: 834
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA 