Scalapack SVDroutine results do not match Lapack svdresults.

[vikas]

Hello,

We have a written a small program to calculate SVD of a matrix using the PDGESVD routine found in ScaLapack 1.8 library. The programs runs and completes the the calculations without any errors.

But when we compare the results obtained by the PDGESVD routine against the results from the dgesvd routine in lapack, we observe a lot of variations between the two results.

Please note there is only one process in the process grid when calling the PDGESVD routine.

It would be of great help if any one could shed some light on what we might be doing wrong in the program. I have pasted the code below for your reference:

Thanks in advance

-Vikas

START OF CODE

CALL BLACS_PINFO( IAM, NPROCS )


*
* Open file to read configs
*
OPEN( NIN, FILE='TestProg1Config.dat', STATUS='OLD' )
READ( NIN, FMT = * ) TOTMEM
READ( NIN, FMT = * ) lwork
READ( NIN, FMT = * ) SIZE
READ( NIN, FMT = * ) NB
READ( NIN, FMT = * ) NPROW
READ( NIN, FMT = * ) NPCOL

MEMSIZ = TOTMEM / CPLXSZ
ALLOCATE (MEM(MEMSIZ), STAT=INFO)

WRITE( *, FMT = * )'Allocation successful:', INFO

M = SIZE
N = SIZE

MYROW = 1
MYCOL = 1


CALL CPU_TIME(STARTTIME)
CALL BLACS_GET( -1, 0, ICTXT )
CALL BLACS_GRIDINIT( ICTXT, 'Row-major', NPROW, NPCOL )
CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )

IA = 1
JA = 1
IU = 1
JU = 1
IVT = 1
JVT = 1

LDA = NUMROC( M, NB, MYROW, 0, NPROW )
LDA = MAX( 1, LDA )
NQ = NUMROC( N, NB, MYCOL, 0, NPCOL )
LDU = LDA
SIZEQ = NUMROC( SIZE, NB, MYCOL, 0, NPCOL )
LDVT = NUMROC( N, NB, MYROW, 0, NPROW )
LDVT = MAX( 1, LDVT )



CALL DESCINIT( DESCA, M, N, NB, NB, 0, 0, ICTXT, LDA, INFO )
CALL DESCINIT( DESCU, M, N, NB, NB, 0, 0, ICTXT, LDU, INFO )
CALL DESCINIT( DESCVT, N, N, NB, NB, 0, 0, ICTXT, LDVT,
$ INFO )




IPA = 1
IPS = IPA + LDA*NQ
IPW = IPS + SIZE
*
IPU = IPW
IPVT = IPW

! Dummy call to get the space needed.
CALL PDGESVD( 'V', 'V', M, N, MEM( IPA ), IA, JA, DESCA,
$ MEM( IPS ), MEM( IPU ), IU, JU, DESCU,
$ MEM( IPVT ), IVT, JVT, DESCVT,
$ MEM( IPW ), -1, DINFO )

WPDGESVD = INT( MEM( IPW ) )

WRITE( *, FMT = * )'Extra Workspace for SVD:', WPDGESVD

IPVT = IPU + LDU*SIZEQ


*
* 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( *, FMT = * ) 'MEMORY'
END IF

CALL PDLAREAD( 'MAT.dat', MEM(IPA), DESCA, 0, 0,MEM(IPW) )

CALL PDGESVD( 'V', 'V', M, N, MEM(IPA), IA, JA, DESCA,
$ MEM(IPS), MEM(IPU), IU, JU, DESCU,
$ MEM(IPVT), IVT, JVT, DESCVT,
$ MEM(IPW), WPDGESVD , INFO )

WRITE( *, FMT = * )
$ 'INFO:',INFO

* Printing the output to file.
CALL PDLAWRITE( 'TestProgSOLU.dat', M, M, MEM(IPU), 1, 1, DESCU,
$ 0, 0, MEM(IPW) )

CALL PDLAWRITE( 'TestProg1SOLVT.dat', M, M,MEM(IPVT), 1, 1,DESCVT,
$ 0, 0, MEM(IPW) )

CALL PDLAWRITE( 'TestProg1SOLS.dat', M, 1, MEM(IPS), 1, 1, DESCU,
$ 0, 0, MEM(IPW) )

CALL BLACS_GRIDEXIT( ICTXT )
CALL BLACS_EXIT( 0 )


END OF CODE

[Julien Langou]

Your code is probably fine. Can you give us a small example of input/output
matrices. Julien.

[vikas]

Hi,

Thanks for your response. I am pasting the input/output that we are observering here.

-Vikas

Input Matrix A(3 X 3)

3 3
1
2
3
4
5
6
7
8
9


Output Matrix:

S (Singular values)

-0.382364417679032578D-01
0.100297609754529837D+01
-0.214837238368396355D+00

U(3 x 3)

3 3
0.113728267941119632D+01
-0.157279033691498593D+00
0.128467061064330867D+00
0.363150454775932219D+00
0.920559634711694885D+00
-0.382364417679032578D-01
0.100297609754529837D+01
-0.214837238368396355D+00
-0.382364417679032578D-01

VT(3 X 3)

3 3
0.100297609754529837D+01
-0.214837238368396355D+00
0.408248290463863850D+00
0.938866323126370861D+00
-0.185567362498326366D+00
0.907839424175862914D-01
-0.213063140037730303D-01
0.422961181732366825D-01
0.807082972555794725D-01

[thombos]

Hey,

I'm a FORTRAN noob, but your code looks to me like the work space, U and S somehow all occupy the same area in memory. This might explain the overlap in your result numbers. You think you could space the memory areas for these things further apart so one can see whats what ?

[vikas]

Nice Catch, overlapping memory could be an issue. Will look in this direction.

Thanks
-Vikas