PDSYEVD gives wrong eigenvalues

[Petros]

Hello!
I have tried to implement the scaLAPACK diagonalization routine PDSYEVD in order to diagonalize a symmetric real matrix (named B in the attached code below).
However, I fail miserably every time I try. The subroutine in which PDSYEVD is used runs without any complaints, however the eigenvalues and eigenvectors always
come out wrong. Below I have listed the code. Here B(N,N) is the real symmetric matrix (global form) that I want to diagonalize. (In the below example the matrix
EIGENVECT is not assigned any output.) The routines GRIDSETUP and BLOCKSET just calculate the computational grid and block size.

For example N = 190 , NPROCS = 24 => NB= 47, NPROW=4, NPCOL=6.

PLEASE SOMEONE HELP ME I AM GOING INSANE FROM TRYING TO MAKE THIS WORK!!




SUBROUTINE diaghscalapack( B,N,EIGENVAL,EIGENVECT)
IMPLICIT NONE
INTEGER :: LWORK, MAXN,LIWORK
DOUBLE PRECISION :: rtmp(4)
INTEGER :: itmp( 4 )
INTEGER :: CONTEXT, I,J, IAM, INFO, MYCOL, MYROW, N, NB, MB, NPCOL, NPROCS, NPROW,LNROWSA,LNCOLSA,SQRTNP
INTEGER :: DESCA(9), DESCZ(9)
DOUBLE PRECISION :: B(N,N),EIGENVAL(N),EIGENVECT(N,N)
DOUBLE PRECISION, ALLOCATABLE :: A(:,:),Z(:,:),WORK(:)
INTEGER, ALLOCATABLE :: IWORK(:)
INTEGER :: NUMROC
EXTERNAL :: BLACS_EXIT, BLACS_GET, BLACS_GRIDEXIT, BLACS_GRIDINFO, BLACS_GRIDINIT, BLACS_PINFO,BLACS_SETUP, DESCINIT, PDLAPRNT,PDSYEV,PDYSEVD,NUMROC,blockset
!
! Initialize the BLACS
!
CALL BLACS_PINFO( IAM, NPROCS )

! Calculating the number of rows, NPROW, and columns, NPCOL of the processor grid:

CALL GRIDSETUP(NPROCS,NPROW,NPCOL)
!
! Initialize a single BLACS context
!
CALL BLACS_GET( -1, 0, CONTEXT )
CALL BLACS_GRIDINIT( CONTEXT, 'R', NPROW, NPCOL )
CALL BLACS_GRIDINFO( CONTEXT, NPROW, NPCOL, MYROW, MYCOL )

! Calculating the block sizes:

CALL BLOCKSET( NB,MB, N, NPROW, NPCOL )

! Calculate the number of rows, LNROWSA, and the number of columns, LNCOLSA, for the local matrices A:

LNROWSA = NUMROC(N,NB,MYROW,0,NPROW)
LNCOLSA = NUMROC(N,NB,MYCOL,0,NPCOL)

ALLOCATE(A(LNROWSA,LNCOLSA),Z(LNROWSA,LNCOLSA))

A = 0.0d0

! Dstribute the global matrix B onto the local matrices A:

CALL DESCINIT( DESCA, N, N, NB, NB, 0, 0, CONTEXT, LNROWSA, INFO )

DO I=1,N
DO J=1,N
IF ( J .GE. I ) CALL PDELSET( A, I, J, DESCA, B(I,J) )
! Since B is supposed to be symetric but only contains
! true matrix elements in the upper triangular part
IF ( J .LT. I ) CALL PDELSET( A, I, J, DESCA, B(J,I) )
ENDDO
ENDDO

CALL DESCINIT( DESCZ, N, N, NB, NB, 0, 0, CONTEXT, LNROWSA, INFO )

! Ask PDSYEVD to compute the entire eigendecomposition

LWORK = -1
LIWORK = 1
rtmp(:) = 0.0d0
itmp(:) = 0

CALL PDSYEVD( 'V', 'U', N, A, 1, 1, DESCA, EIGENVAL, Z, 1, 1, DESCZ, rtmp, LWORK,itmp, LIWORK, INFO )

LWORK = MAX( 131072, 2*INT( rtmp(1) ) + 1 )
LIWORK = MAX( 8*N , itmp(1) + 1 )
ALLOCATE( IWORK(LIWORK),WORK(LWORK))

CALL PDSYEVD( 'V', 'U', N, A, 1, 1, DESCA, EIGENVAL, Z, 1, 1, DESCZ, WORK, LWORK,IWORK, LIWORK, INFO )

! Print out the eigenvectors

CALL PDLAPRNT( N, N, Z, 1, 1, DESCZ, 0, 0, 'Z', 6, WORK )

CALL BLACS_GRIDEXIT( CONTEXT )

CALL BLACS_EXIT( 1 )

END SUBROUTINE diaghscalapack

[Petros]

Hello!
I found the error my self! The error was pretty basic/stupid and was due to that the NxN symmetric matrix B was not accessible to
all the processors (ranks) prior to the call of the subroutine diaghscalapack, it existed only on rank=0. By making sure that B existed on
the entire grid the problem was solved. Perhaps not the best solution from the perspective of memory use, since only limited local parts need
to be accessible on each of the processors, but hey a quick fix is also a fix.

Hope that this can be of some help to anyone of you brave ScaLapack users out there!

[Julien Langou]

Hi Petros, thanks for sharing your experience with ScaLAPACK! Julien.