PZGESVD routine

[beginner]

Dear all I am using scalapack to perform a singular value decomposition on a big matrix

I compare the results from ScaLAPACK with the LAPACK code that I have and I see that when I use the PZGESVD routine for singular value decomposition the
routine gives the same singular values as the LAPACK routine ZGESVD but I obtain different left and right singular vectors U and VT.

Does any body experience a similar problem?
Here is the piece of code that I wrote

ESIZE = MIN( M, N )
allocate( S(ESIZE) )

LDA = NUMROC( M , MB , MYROW , 0 , NPROW )
LDA = max( 1, LDA )

LDU = LDA

LDVT = NUMROC( N , MB , MYROW , 0 , NPROW )
LDVT = max( 1, LDVT )

CALL DESCINIT( DESCA , M , N , MB , NB ,
& RSRC , CSRC , ICTXT , LDA , INFO )
CALL DESCINIT( DESCU , M , M , MB , NB ,
& RSRC , CSRC , ICTXT , LDU , INFO )
CALL DESCINIT( DESCVT , N , N , MB , NB ,
& RSRC , CSRC , ICTXT , LDVT , INFO )

A = M/NPROW + MB
B = N/NPROW + NB
allocate( A1(A,B) , A2(A,B) )
allocate( WORKA1( A + MB + 1 ) , WORKA2( A + MB + 1 ) )

A = M/NPROW + MB
B = M/NPROW + NB
allocate( U(A,B) )
allocate( WORKU( A + MB + 1 ) )

A = N/NPROW + MB
B = N/NPROW + NB
allocate( VT(A,B) )
allocate( WORKVT( A + MB + 1 ) )

CALL DASSIGN_MATRIX( A1 , A2 , M , N , DESCA ,
& 0 , 0 , WORKA1 , WORKA2 )

LWORK1 = -1
ALLOCATE( WORK(1) , RWORK(1) )
CALL PZGESVD( 'V' , 'V' , M , N , A1 , 1 , 1 , DESCA , S ,
& U , 1 , 1 , DESCU ,
& VT , 1 , 1 , DESCVT,
& WORK , LWORK1 , RWORK , INFO )

LWORK1 = INT(dble(WORK(1)))
LWORK2 = INT(dble(RWORK(1)))

DEALLOCATE( WORK , RWORK )

ALLOCATE( WORK(LWORK1) , RWORK(LWORK2) )
CALL PZGESVD( 'V' , 'V' , M , N , A1 , 1 , 1 , DESCA , S ,
& U , 1 , 1 , DESCU ,
& VT , 1 , 1 , DESCVT,
& WORK , LWORK1 , RWORK , INFO )
DEALLOCATE(WORK,RWORK)

[admin]

Hi,

1) for two distinct singular values, singular vector does not have to be the same up to multiplication with a scalar number of the kind e^(i.theta). While this is a simple +/- in the real numbers, in the complex numbers it is hard to see that to vectors are multiple by e^(i.theta). So please check this.

2) if this is singular values are similar (or very close), then the
subspace is invariant, it's even harder to compare the answer of LAPACK and ScaLAPACK.

3) All in all, you can check that
a) || A - U * S * V^H || / || A || is small
b) || I - U * U^H || is small
c) || I - V * V^H || is small
for both LAPACK and ScLAPACK. They both have A correct answer. There is no such thing has THE correct answer for this problem.

Hope it helps
Julie

[beginner]

Thanks for the reply.

I will try the checks at first.

[beginner]

I have been looking in the results and I could see that they are depending on the value that I choose for MB and NB meaning the block sizes.
Does anyone experience also that?

[Julien Langou]

From the algorithm used, I think the results should be "identical" for two different pairs of block sizes ...
Can you try the check please? || A - U * S * V^H || / || A || and orthogonality of U and V?
--JL

[beginner]

I did the test for the PZGESVD output and the results is that the output matrix is not orthogonal.

Here is the code I used for the test

! --- Generate S matrix

LDSM = NUMROC( ESIZE , MB , MYROW , 0 , NPROW )
CALL DESCINIT( DESCSM , ESIZE , ESIZE , MB , NB ,
& RSRC , CSRC , ICTXT , LDSM , INFO )
A = LDSM
B = NUMROC( ESIZE , NB , MYCOL , 0 , NPCOL )
allocate( SM(A,B) )
allocate( WORKSM( MB ) )

CALL DVECTOR2CMATRIXN( S , SM , ESIZE , ESIZE , DESCSM ,
& 0 , 0 , WORKSM )

! --- Singular value decomposition check

LDC1 = NUMROC( M , MB , MYROW , 0 , NPROW )
CALL DESCINIT( DESCC1 , M , ESIZE , MB , NB ,
& RSRC , CSRC , ICTXT , LDC1 , INFO )
A = LDC1
B = NUMROC( ESIZE , NB , MYCOL , 0 , NPCOL )
allocate( C1(A,B) )

LDC = NUMROC( M , MB , MYROW , 0 , NPROW )
CALL DESCINIT( DESCC , M , N , MB , NB ,
& RSRC , CSRC , ICTXT , LDC , INFO )
A = LDC
B = NUMROC( N , NB , MYCOL , 0 , NPCOL )
allocate( C(A,B) )

CALL PZGEMM( 'N' , 'N' , M , ESIZE , ESIZE , ALPHA ,
& U , 1 , 1 , DESCU ,
& SM , 1 , 1 , DESCSM ,
& BETA ,
& C1 , 1 , 1 , DESCC1 )

CALL PZGEMM( 'N' , 'N' , M , N , ESIZE , ALPHA ,
& C1 , 1 , 1 , DESCC1 ,
& VT , 1 , 1 , DESCVT ,
& BETA ,
& C , 1 , 1 , DESCC )

CALL IDENTITY_MATRIX( I , M , N , DESCA , 0 , 0 , WORKA1 )

CALL PZGEMM( 'N' , 'N' , M , N , ESIZE , ALPHA1 ,
& A1C , 1 , 1 , DESCA ,
& I , 1 , 1 , DESCA ,
& BETA1 ,
& C , 1 , 1 , DESCC )

iroffc = mod(0,MB)
mp0 = NUMROC( M + iroffc , MB , MYROW , 0 , NPCOL )
allocate( WORKC(mp0) )

NORM1 = PZLANGE( 'I' , M , N , C , 1 , 1 , DESCC , WORKC )
NORM2 = PZLANGE( 'I' , M , N , A1C , 1 , 1 , DESCA , WORKC )
RESID = NORM1 / NORM2

! --- Check ortogonality U ( M ESIZE ) VT ( ESIZE N )

LDUC = NUMROC( M , MB , MYROW , 0 , NPROW )
CALL DESCINIT( DESCUC , M , M , MB , NB ,
& RSRC , CSRC , ICTXT , LDUC , INFO )
A = LDUC
B = NUMROC( M , NB , MYCOL , 0 , NPCOL )
allocate( UC(A,B) )

iroffuc = mod(0,MB)
mp0uc = NUMROC( M + iroffc , MB , MYROW , 0 , NPCOL )
allocate( WORKUC(mp0uc) )

LDVTC = NUMROC( N , MB , MYROW , 0 , NPROW )
CALL DESCINIT( DESCVTC , N , N , MB , NB ,
& RSRC , CSRC , ICTXT , LDVTC , INFO )
A = LDVTC
B = NUMROC( N , NB , MYCOL , 0 , NPCOL )
allocate( VTC(A,B) )

iroffvtc = mod(0,MB)
mp0vtc = NUMROC( ESIZE + iroffc , MB , MYROW , 0 , NPCOL )
allocate( WORKVTC(mp0vtc) )

CALL PZGEMM( 'N' , 'T' , M , M , ESIZE , ALPHA ,
& U , 1 , 1 , DESCU ,
& U , 1 , 1 , DESCU ,
& BETA ,
& UC , 1 , 1 , DESCUC )

NORM3 = PZLANGE( 'I' , M , M , UC , 1 , 1 , DESCU , WORKUC )

CALL PZGEMM( 'T' , 'N' , N , N , ESIZE , ALPHA ,
& VT , 1 , 1 , DESCVT ,
& VT , 1 , 1 , DESCVT ,
& BETA ,
& VTC , 1 , 1 , DESCVTC )

NORM4 = PZLANGE( 'I' , N , N , VTC , 1 , 1 ,
& DESCVTC , WORKVTC )

WRITE(*,'(a,1x,I0,3(ES15.6,1x))')
&'Residual', IAM, RESID , NORM3 , NORM4

The output is the following

Residual 0 2.447124E-08 2.780209E+01 1.000000E+00
Residual 1 2.447124E-08 2.780209E+01 1.000000E+00
Residual 2 2.447124E-08 2.780209E+01 1.000000E+00
Residual 3 2.447124E-08 2.780209E+01 1.000000E+00

[Julien Langou]

Well you did the hardest ... the RESID = || A - USVT || / || A || checks look OK, correct?

For your orthogonality checks:
1) when do you initialize UC to identity? Don't you want UC = Identity?
2) In the line: CALL PZGEMM( 'N' , 'T' , M , M , ESIZE , ALPHA ,
you want the second 'T' to be a 'C' !!!

I remind you: you want to check:
|| I - U^H * U ||
( or || I - U * U^H || , this is the same when the matrix U is square)

Julien.

[beginner]

Thanks for the help once again, i wrote the part for the orthogonality once again and this are the results I get on 8 processors
The result does not change when I change the number of processors and the block size.

IAM = 0 SVD Residual = 2.447123E-08
IAM = 2 SVD Residual = 2.447123E-08
IAM = 1 SVD Residual = 2.447123E-08
IAM = 6 SVD Residual = 2.447123E-08
IAM = 5 SVD Residual = 2.447123E-08
IAM = 3 SVD Residual = 2.447123E-08
IAM = 7 SVD Residual = 2.447123E-08
IAM = 4 SVD Residual = 2.447123E-08
IAM = 0 U check = 1.791979E-13
IAM = 1 U check = 1.791979E-13
IAM = 2 U check = 1.791979E-13
IAM = 6 U check = 1.791979E-13
IAM = 5 U check = 1.791979E-13
IAM = 3 U check = 1.791979E-13
IAM = 7 U check = 1.791979E-13
IAM = 4 U check = 1.791979E-13
IAM = 0 VT check = 1.409582E-14
IAM = 2 VT check = 1.409582E-14
IAM = 5 VT check = 1.409582E-14
IAM = 6 VT check = 1.409582E-14
IAM = 1 VT check = 1.409582E-14
IAM = 3 VT check = 1.409582E-14
IAM = 7 VT check = 1.409582E-14
IAM = 4 VT check = 1.409582E-14

I can assume now that my Singular Value Decomposition should be right. I compared anyway the results with my serial program written
using LAPACK and the singular right and left vectors are still different, maybe they are normalized in a different way ?

In my next step for the algorithm I need to compute eigenvalues and eigenvectors of a non symmetric matrix.
I had a look at the subroutine but the parallel version of the subroutine zgeev seems not to be available anymore.

[Julien Langou]

I can assume now that my Singular Value Decomposition should be right.

Yes it is. Congratulations. (For the successful use of the subroutines but also for writing the check correctly.)

I compared anyway the results with my serial program written
using LAPACK and the singular right and left vectors are still different, maybe they are normalized in a different way ?

Please reread our initial answer. There is not unicity of the singular vectors.

• In exact arithmetic, if the singular value is unique, then singular vectors may differ by a e^(i.theta) (so a complex scalar of modulus 1). In real arithmetic, this is a sign, so that's easy to spot. In complex arithmetic, this is really harder to check with your eyes that two vectors differ by e^(i.theta). What is unique is the invariant subspace associated with this singular value: so any vector of norm 1 in that subspace of dimension 1 is fine. Although the singular vectors are different, there is no wrong answer. LAPACK is correct, ScaLAPACK is correct.
  
• In exact arithmetic, if the singular value is multiple say of multiplicity 4 (e.g.): what is unique here is once more the invariant subspace of dimension 4 associated with this singular value: so any orthogonal basis of this subspace of dimension 4 is fine. This gives much more freedom to the singular value solvers ... Once more, although the singular vectors are different, there is no wrong answer. LAPACK is correct, ScaLAPACK is correct.
  
• In 32-bit (or 64-bit) arithmetic, then if the singular value are close together ... well you end up with the same problem than in exact arithmetic with a multiple singular value. I think this makes sense but it is hard to describe in this post what's going on precisely.

In my next step for the algorithm I need to compute eigenvalues and eigenvectors of a non symmetric matrix.
I had a look at the subroutine but the parallel version of the subroutine zgeev seems not to be available anymore.

Oups ... This is coming, this is coming ... maybe register yourself at lapack-announce@cs.utk.edu. We hope to have the feature available within the next 6 months. I do not recommend you to try anything in ScaLAPACK at this point. Sorry for the bad news.

Cheers,
Julien.

[beginner]

Thanks once again for the reply.
I found on the internet some subroutines that should do the work for the eigenvalue eigenvectors problem but I do not know whether they were deeply tested.
On the other side since the SVD reduces the size of the matrix I have to work on i simply use the lapack subroutine ZGEEV on a single processors and then i send the results to the others.

It would be nice anyway to have also the eigenvalue part in a parallel version so i will be looking forward for the update.