Page 1 of 1

bug with dgesdd

PostPosted: Sat Jun 19, 2010 10:39 am
by Bruno

in the octave bug list there is a matrix example where dgesdd performs badly
(see ... 87#options).
I don't know if lapack team is aware of this problem (at least I didn't see it
on this forum).

I have transformed the matrix in a ascii format and written a small f90 program
which reads it and use dgesdd or dgesvd but I am unable to attach
the 2 files (what is possible to attach to a message for this forum ?) So I
have put the 2 files in my homepage :

On the given matrix :

1/ dgesvd works quite well, I got :

|| U*diag(sigma)*VT - A ||_1 / || A ||_1 near 10^(-14)

2/ but with dgesdd, first the last singular value is negative !!! and
the previous relative error is about 48.4 !

I have tested using 2 differents blas lib and got the same results.


Re: bug with dgesdd

PostPosted: Mon Jun 21, 2010 12:10 pm
by admin
Thanks Bruno!
We actually knew that there was a serious bug in dgesdd but we couldn't find a matrix to reproduce it.
For the current problems found in LAPACK, you can always refer to
The dgesdd problem is bug0025 as you can see.
The good news is that I was able to reproduce it and so I confirm the problem.
Your code and your matrix are really great and straightforward, so we will start debugging the dgesdd routine.
We keep you posted.
Merci beaucoup
Julie et Julien Langou

Re: bug with dgesdd

PostPosted: Mon Jun 21, 2010 2:46 pm
by Bruno
The matrix has been posted on the octave bug tracker
by Jaroslav Hajek, I have only translated it in ascii and
wrote the small f90 code.

Thanks for the ref bug,


Re: bug with dgesdd

PostPosted: Mon Jun 21, 2010 4:56 pm
by Julien Langou

Thanks for the bug report!

So I have taken a look at all this.

Good news: I think I have found where the bug is.
Bad news: I think the routine should have failed with INFO = 1.
So does not really help. But that helps somewhat!

So line 353 from 374 of DBDSC reads:
Code: Select all
            IF( ICOMPQ.EQ.2 ) THEN
               CALL DLASD0( NSIZE, SQRE, D( START ), E( START ),
     $                      U( START, START ), LDU, VT( START, START ),
     $                      LDVT, SMLSIZ, IWORK, WORK( WSTART ), INFO )
     $                      E( START ), Q( START+( IU+QSTART-2 )*N ), N,
     $                      Q( START+( IVT+QSTART-2 )*N ),
     $                      IQ( START+K*N ), Q( START+( DIFL+QSTART-2 )*
     $                      N ), Q( START+( DIFR+QSTART-2 )*N ),
     $                      Q( START+( Z+QSTART-2 )*N ),
     $                      Q( START+( POLES+QSTART-2 )*N ),
     $                      IQ( START+GIVPTR*N ), IQ( START+GIVCOL*N ),
     $                      N, IQ( START+PERM*N ),
     $                      Q( START+( GIVNUM+QSTART-2 )*N ),
     $                      Q( START+( IC+QSTART-2 )*N ),
     $                      Q( START+( IS+QSTART-2 )*N ),
     $                      WORK( WSTART ), IWORK, INFO )
               IF( INFO.NE.0 ) THEN
               END IF
            END IF

and I strongly believe the check on INFO should be outside the ELSE statement.

In your case, DBDSDC is calling DLASD0 twice (well with vecLib not with reference BLAS !!!), the first call returns an INFO = 1 so everything should stop. The second call returns INFO = 0 and overrides the first INFO. This is wrong.

So I'll send this to the lapackers and hopefully we will patch this soon.


Re: bug with dgesdd

PostPosted: Fri Jun 25, 2010 8:17 am
by Bruno
Thanks Julien. I have another question : when dgesdd "normally" outputs
with info=0, can we considered that the singular values and vectors have
near about the same accuray than with dgesvd ? (is there some paper
about the stability of dgessd ?)

Cordiales salutations

Re: bug with dgesdd

PostPosted: Thu Jul 01, 2010 2:57 am
by Julien Langou
Hello, update: Last week, I forwarded your question around because I do not know the answer for sure. I did not receive any answer yet. Bottom line there should be not much difference in term of the quality of the answer. The first step for both is a reduction to bidiagonal form. The second step whether you use GESDD or GESVD has about the same accuracy or more accurate than the first. So at the end, you'll get the same accuracy more or less. Cheers, Julien.

Re: bug with dgesdd

PostPosted: Mon Jul 12, 2010 3:56 am
by Bruno
Thanks for the information Juien (and sorry for the delay I suppose
I forgot the dgesdd pb...).

Finally this means that for a singular elements computations one can
try first dgesdd and then dgesvd in case the first fails.


Confused with square root of matrix and dgesdd results

PostPosted: Tue Aug 03, 2010 7:19 pm
by vukicevic
Hi. I'm trying to find square root of the folowing matrix:

A =

I use dgesdd function, and I've got:

U =
(-0.769509, 0.638636)
(0.638636, 0.769509)

S =
(22.2788, 0)
(0, 10.2788)

VT =
(-0.769509, 0.638636)
(-0.638636, -0.769509)

Until now everything seems ok, because

U * S * VT = A =

Now it should be possible to get square root of A simply by taking square roots of eigenvalues in matrix S:

S**(1/2) =
(4.72004, 0)
(0, 3.20606)

and square root of A should be:

A**(1/2) = U * S**(1/2) * VT

But when I had tried to check the result, I've got:

A**(1/2) * A**(1/2) =
(17.3846, -5.89723)
(-5.89723, 15.1731)

which is obviously wrong.

I'm a litle confused with that result because the same procedure works fine with e.g.

B =

I suppose that I can try to solve the problem with dsyevr function, but it is interesting to know what is wrong.

Best regards, Vlatko

Re: bug with dgesdd

PostPosted: Wed Aug 04, 2010 2:35 pm
by Julien Langou
You will not be able to compute the square root of a matrix using an SVD decomposition!