LAPACK Archives

[Lapack] Fwd: Re: incorrect singular values in DGESDD ?

Sorry, I forgot to cc the group in my response to the user.
Osni


-------- Forwarded Message --------
Subject:        Re: [Lapack] incorrect singular values in DGESDD ?
Date:   Mon, 14 Sep 2015 14:15:02 -0700
From:   Osni Marques <oamarques@Domain.Removed>
To:     Bernhard Dick <Bernhard.Dick@Domain.Removed>



Dear Bernhard,

It looks like there is a typo in the test program TestSVD.for, lines 68-79:

C
C   do a work space query
C
       LWORK = -1
       CALL DGESDD ( CJOB,  M, N, A, MAXN, S, U, MAXN, VT, MAXN,
      $                WORK, LWMAX, IWORK, INFO )
       LWORK = INT(WORK(1)+0.5D0)
       WRITE (*,*) ' opt workspace',LWORK,LWMAX
       IF (LWORK.GT.LWMAX .OR. IWDIM.GT.IWMAX) STOP
       CALL DGESDD ( CJOB,  M, N, A, MAXN, S, U, MAXN, VT, MAXN,
      $                WORK, LWMAX, IWORK, INFO )
       WRITE (*,*) 'Testing DGESDD, INFO=',INFO

The first call to DGESDD should be with LWORK in the 12th argument 
(instead of LWMAX), otherwise DGESDD does perform the SVD of A. Upon 
return the values in A are destroyed. In addition, with CJOB='0', U is 
overwritten on A. Given the input matrix (i.e. it is diagonal), U is the 
identity matrix. This matrix  is then passed to DGESDD in the second 
call, which returns singular values equal to one (as expected).

Best regards,

Osni


On 9/12/2015 7:34 AM, Bernhard Dick wrote:
Dear Sir or Madam,
I am using Lapack for more than 10 years, but now I get quite confused 
by the observation that I do not get correct singular values from 
DGESDD with the JOBZ='O' option. This troubles me in particular since 
I was using this subroutine quite a lot in previous programs, and they 
seemed to work. I should say that I am using the OpenWatcom compiler 
on Windows, so I built my own libraries for Lapack and BLAS, the 
latest versions were Lapack 3.1.1 and now Lapack3.5.0. So I was 
suspecting that something is wrong with the Watcom compiler. But now I 
used the ABsoft FORTRAN compiler which uses its own built-in Lapack, 
with the same result.
I attach a test program that shows the problem: A test matrix is set 
up which is zero except for the diagonal, which has ascending numbers 
A(I,I) = DFLOAT(I). So these should also be the singular values, 
and the singular vectors are all of the form (0..0,1,0..0). I compare 
DGESDD with DGESVD. Input is NROW= number of rows in A, NCOL=number of 
columns in A, JOBZ=either 'A' or 'O' (the quotes are needed in the 
input with the WATCOM compiler, but not with the ABsoft compiler). I 
observe the following:
DGESVD gives the correct result for all NROW >= NCOL
DGESDD gives the same correct results only for JOBZ='A'
For JOBZ='O', however, all singular values produced by DGESDD are 1.00000
The same is observed with ABsoft FORTRAN and built in Lapack, and 
Watcom-FORTRAN and my only build of Lapack and BLAS. So I would 
conclude that the problem is not with the compilers, unless they do 
both the same wrong thing. I disabled optimization in the WATCOM case.
I include the source and the two executables. I hope you see something 
that I am doing wrong, because otherwise I am completely lost. I am 
particularly interested in the 'O' option because in this case A and U 
can be physically the same memory location.
with best regards,
Bernhard Dick
p.s.:  the executables were rejected by the mail server. If you would 
like to have them I can send them encrypted.
Prof. Dr. Bernhard Dick
Institut f?r Physikalische und Theoretische Chemie
Universit?t Regensburg
D 93040 Regensburg
Tel.: +49 941 943 4487
Fax: +49 941 943 4488
www-dick.chemie.uni-regensburg.de


_______________________________________________
Lapack mailing list
Lapack@Domain.Removed
http://lists.eecs.utk.edu/mailman/listinfo/lapack



-------------- next part --------------
An HTML attachment was scrubbed...
URL: 
<http://lists.eecs.utk.edu/mailman/private/lapack/attachments/20150914/72d7c405/attachment.html>

<Prev in Thread] Current Thread [Next in Thread>
  • [Lapack] Fwd: Re: incorrect singular values in DGESDD ?, oamarques <=


For additional information you may use the LAPACK/ScaLAPACK Forum.
Or one of the mailing lists, or