EISPACK vs. SCALAPACK different results..??

[2linear]

Hi There,
I am working in a code that computes the eigenvalues and eigenvectors using PSSYEV routine. The code has been tested/verified using the different matrices provided in the Scalapack examples website and the results (eigenvalues and eigenvectors) are exactly the same reported on the website. My ultimate goal is to use the code to solve the eigenproblem for large matrices. Right now I am working with a matrix of dimension 4455 x 4455 and the code runs fast (a couple of minutes with 4 CPUs). However, the results I am getting are way different compared to those obtained using EISPACK. I have used EISPACK in the past so there is not problem with the EISPACK code. I would like to know if It is possible that PSSYEV doesn't handle properly such large matrices, or, if there is something wrong with the code itself. Any help on this is much appreciated.
The code
===========================================================================

Code: Select all

      PROGRAM TEST
       INTEGER :: irec,i,j,l,k,LWORK,MAXN,LDA
      REAL :: ZERO =0.,MONE=-1.
      INTEGER MEMSIZ, TOTMEM
      INTEGER, PARAMETER :: REALSZ = 4
      INTEGER, PARAMETER :: INTGSZ = 4
      INTEGER, PARAMETER :: BLOCK_CYCLIC_2D = 1
      INTEGER, PARAMETER :: DLEN_ = 9, DT_ = 1,CTXT_ = 2, M_ = 3
      INTEGER, PARAMETER :: N_ = 4, MB_ = 5, NB_ = 6, RSRC_ = 7
      INTEGER, PARAMETER :: CSRC_ = 8, LLD_ = 9
      REAL, PARAMETER :: ONE = 1.0E+0
      CHARACTER*80 ::   OUTFILE
      INTEGER :: IAM, ICTXT, INFO, IPA, IPACPY, IPB, IPPIV, IPX
      INTEGER :: IPW, LIPIV, MYCOL, MYROW, N, NB, NOUT, NPCOL
      INTEGER :: NPROCS, NPROW, NP, NQ, NQRHS, NRHS, WORKSIZ
      REAL :: ANORM, BNORM, EPS, XNORM, RESID
      REAL :: slambda
      INTEGER :: DESCA( DLEN_ ), DESCB( DLEN_ ), DESCX( DLEN_ )
      INTEGER :: DESCZ(DLEN_)
      REAL, DIMENSION(:), ALLOCATABLE :: MEM
      REAL, DIMENSION(:,:), ALLOCATABLE :: Z
      REAL, DIMENSION(:), ALLOCATABLE :: W
      REAL, DIMENSION(:), ALLOCATABLE :: WORK

      EXTERNAL :: BLACS_EXIT, BLACS_GET, BLACS_GRIDEXIT
      EXTERNAL :: BLACS_GRIDINFO, BLACS_GRIDINIT
      EXTERNAL :: BLACS_SETUP, PDLAMODHILB, PSSYEV
      EXTERNAL :: BLACS_PINFO, DESCINIT, IGSUM2D, PDSCAEXINFO, PSGESV
      EXTERNAL :: PSGEMM, PSLACPY, PSLAPRNT, PSLAREAD, PSLAWRITE
*     ..
*     .. External Functions ..
      INTEGER :: ICEIL, NUMROC
      REAL :: PSLAMCH, PSLANGE
      EXTERNAL :: ICEIL, NUMROC
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC :: DBLE, MAX, ABS
        LWORK=4455*4455*8
        MAXN=4455
        LDA=MAXN
        ALLOCATE(WORK(LWORK))
        ALLOCATE(Z(LDA,LDA))
        ALLOCATE(W(LDA))

        TOTMEM = 495437044
*        TOTMEM = 2000000
        MEMSIZ = TOTMEM / REALSZ
         ALLOCATE(MEM(MEMSIZ))
*     Get starting information
*
      CALL BLACS_PINFO( IAM, NPROCS )
      CALL PDSCAEXINFO( OUTFILE, NOUT, N, NRHS, NB, NPROW, NPCOL, MEM,
     $                  IAM, NPROCS )
*
*     Define process grid
*
      CALL BLACS_GET( -1, 0, ICTXT )
      CALL BLACS_GRIDINIT( ICTXT, 'Row-major', NPROW, NPCOL )
      CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
*
*     Go to bottom of process grid loop if this case doesn't use my
*     process
*
      IF( MYROW.GE.NPROW .OR. MYCOL.GE.NPCOL )
     $   GO TO 20
*
      NP    = NUMROC( N, NB, MYROW, 0, NPROW )
      NQ    = NUMROC( N, NB, MYCOL, 0, NPCOL )
*
*     Initialize the array descriptor for the matrix A and B
*
      CALL DESCINIT( DESCA, N, N, NB, NB, 0, 0, ICTXT, MAX( 1, NP ),
     $               INFO )
      CALL DESCINIT( DESCZ, N, N, NB, NB, 0, 0, ICTXT, MAX( 1, NP ),
     $               INFO )


      IF( IAM.EQ.0 ) THEN
       PRINT*,'....1.....THIS..PART..IS.OK'
      ENDIF
*
*     Assign pointers into MEM for SCALAPACK arrays, A is
*     allocated starting at position MEM( 1 )
*
      IPA = 1
      IPACPY = IPA + DESCA( LLD_ )*NQ
      IPB = IPACPY + DESCA( LLD_ )*NQ
      LIPIV = ICEIL( INTGSZ*( NP+NB ), REALSZ )
      IPW = IPB+MAX( NP, LIPIV )
*
      WORKSIZ = NB
*
*     Check for adequate memory for problem size
*
      INFO = 0
      IF( IPW+WORKSIZ.GT.MEMSIZ ) THEN
         IF( IAM.EQ.0 )
     $      WRITE( NOUT, FMT = 9998 ) 'test', ( IPW+WORKSIZ )*REALSZ
         INFO = 1
      END IF
*
*     Check all processes for an error
*
      CALL IGSUM2D( ICTXT, 'All', ' ', 1, 1, INFO, 1, -1, 0 )
      IF( INFO.GT.0 ) THEN
         IF( IAM.EQ.0 )
     $      WRITE( NOUT, FMT = 9999 ) 'MEMORY'
         GO TO 10
      END IF
*
*     Read from file and distribute matrix A
*
      CALL PSLAREAD( 'test4455x4455.dat', MEM( IPA ), DESCA, 0, 0,
     $               MEM( IPW ) )
*
*     Make a copy of A for checking purposes
*
      IF( IAM.EQ.0 ) THEN
         PRINT*,'....2..THIS...PART...OK..'
      ENDIF

      CALL PSLACPY( 'All', N, N, MEM( IPA ), 1, 1, DESCA,
     $              MEM( IPACPY ), 1, 1, DESCA )
      IF( IAM.EQ.0 ) THEN
         PRINT*,'....3..THIS...PART...OK..'
       ENDIF

      CALL PSSYEV( 'V', 'U', N, MEM(IPA), 1, 1, DESCA, W, Z, 1, 1,
     $             DESCZ, MEM(IPW), LWORK, INFO )


      CALL PSLAPRNT( N, N, Z, 1, 1, DESCZ, 0, 0, 'Z', 6, MEM(IPW) )
      CALL PSLAWRITE( 'EIGVSOL.dat', N, N, Z, 1, 1, DESCZ,
     $                0, 0, MEM( IPW ) )

      IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
         PRINT *, ' N =', N
                slambda = 0.0
                do 60 l=1,N
                   slambda = slambda + W(l)
  60            continue
         DO 100 I = N, N-10,-1
            PRINT *, ' W(', I, ')=', W(I)/slambda*100., ';'
 100    CONTINUE
      END IF

   10 CONTINUE
*
      CALL BLACS_GRIDEXIT( ICTXT )
*
   20 CONTINUE
*
*     Print ending messages and close output file
*
      IF( IAM.EQ.0 ) THEN
         WRITE( NOUT, FMT = * )
         WRITE( NOUT, FMT = * )
         WRITE( NOUT, FMT = 9997 )
         WRITE( NOUT, FMT = * )
         IF( NOUT.NE.6 .AND. NOUT.NE.0 )
     $      CLOSE ( NOUT )
      END IF
*
      CALL BLACS_EXIT( 0 )
*
 9999 FORMAT( 'Bad ', A6, ' parameters: going on to next test case.' )
 9998 FORMAT( 'Unable to perform ', A, ': need TOTMEM of at least',
     $        I11 )
 9997 FORMAT( 'END OF TESTS.' )
*
      STOP
*
*     End of TEST
*
      END


[Julien Langou]

Your code looks alright. Why don't you compute the eigenvalues with matlab to check which is correct? J.

[2linear]

Hi Julien,
Thanks a lot for your response. I will try with MATLAB and see what happens.

[2linear]

Julien,

Before moving to the large matrix example I decided to compare the results (for small matrices) between matlab and scalapack to be safe. I am getting strange results, however. I really don't know what went wrong;

1) If I run the scalapack example sample_pssyev_call.f in the website I got the following:

armngrd@arxt27:~/local2/SCALAPACK_INSTALLER/scalapack_installer_0.96/SCALAPACK_INSTALACION/build/scalapack-1.8.0/EXAMPLE$ mpirun -np 4 sample_pssyev_call

Code: Select all

Z(     1,     1)=     -0.485377931880384417D-01
Z(     2,     1)=     -0.122222718667358718D+00
Z(     3,     1)=     -0.282485135303397472D+00
Z(     4,     1)=      0.950214627337748641D+00
Z(     1,     2)=      0.912072227031435201D-01
Z(     2,     2)=      0.426620092092789727D+00
Z(     3,     2)=     -0.877038303257523966D+00
Z(     4,     2)=     -0.201197301593937228D+00
Z(     1,     3)=      0.495193043045525583D+00
Z(     2,     3)=     -0.798642041261893532D+00
Z(     3,     3)=     -0.298844131980119054D+00
Z(     4,     3)=     -0.166273644422072181D+00
Z(     1,     4)=      0.862617629821305298D+00
Z(     2,     4)=      0.406482218545086538D+00
Z(     3,     4)=      0.248391118466086525D+00
Z(     4,     4)=      0.170190725350426925D+00
  N =           4
  A =  hilb(N) + diag([1:-1/N:1/N])
  W(           1 )=  0.30481403606055119      ;
  W(           2 )=  0.58196120214201874      ;
  W(           3 )=  0.90849811000512559      ;
  W(           4 )=   2.3809171279827805      ;
  backerror = A - Z * diag(W) * Z'
  resid = A * Z - Z * diag(W)
  ortho = Z' * Z - eye(N)
  norm(backerror)
  norm(resid)
  norm(ortho)

2) When I run matlab I obtain the following;

Code: Select all

>> [V,D]=eig(hilb(4)+diag([0.25:-0.25:0.25]))

V =

    0.0337    0.2111   -0.6984    0.6830
   -0.3551   -0.7667    0.2255    0.4850
    0.7992    0.0174    0.4419    0.4071
   -0.4839    0.6061    0.5158    0.3640


D =

    0.0001         0         0         0
         0    0.0097         0         0
         0         0    0.2695         0
         0         0         0    2.3968



As you can see the results are different. Do you obtain the same inconsistencies ?. Is there any post-processing applied on the eigenvalues/eigenvectors (e.g. weights) in Matlab ? Or am I doing something nasty ?
Any help on this is much appreciated.
Thanks

[2linear]

Julien,
I realized I made a mistake in the input matrix. It should be hilb(4)+diag([1:..... instead of hilb(4)+diag([0.25...
This machine has not installed matlab so I cannot repeat the experiment with the modified matrix. I will post later with the results.

[2linear]

Hi,
I already tested the modified matrix and the results with MATLAB are exactly the same of those obtained with Scalapack. Now i can move to the large matrix case. I will post later

[2linear]

Hi,
This problem has been fixed. I was having different results between SCALAPACK and EISPACK (for the large matrix case) because there were some "subtle" differences between the input matrix I was using to drive each code. My bad :(. I fixed this problem and now everything is OK. When I was doing the tests for the small matrices I was using the same input matrices to drive the two codes. For that reason I was not having problem with the small matrix case.