LAPACK trunk; issue with ?GESVD

Open discussion regarding features, bugs, issues, vendors, etc.

LAPACK trunk; issue with ?GESVD

Postby sguazt » Fri Nov 18, 2011 7:22 am

Hello,

Just checked out the last SVN version of LAPACK.

The issue I have is related to the computation of workspace size in ?GESVD subroutines

Here below is a test case ZGESVD (file: test_zgesvd.f):

Code: Select all
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        (NIN=5,NOUT=6)
      INTEGER          MMAX, NMAX
      PARAMETER        (MMAX=10,NMAX=24)
      INTEGER          LDA, LDVT
      PARAMETER        (LDA=MMAX,LDVT=NMAX)
*     .. Local Scalars ..
      INTEGER          INFO,  M, N
*     .. Local Arrays ..
      COMPLEX *16      A(LDA,NMAX), DUMMY(1,1), WORK(1,1)
      DOUBLE PRECISION RWORK(5*MIN(MMAX,NMAX))
*     .. Executable Statements ..
      WRITE (NOUT,*) 'F08KPF Example Program Results'
      WRITE (NOUT,*)
*     Skip heading in data file
      READ (NIN,*)
      READ (NIN,*) M, N
      WRITE (NOUT,*) 'READ>> M: ', M, ' - N: ', N
      IF (M.LE.MMAX .AND. N.LE.NMAX) THEN
*
*        Read the m by n matrix A from data file
*
         READ (NIN,*) ((A(I,J),J=1,N),I=1,M)
*
*        Compute the singular values and left and right singular vectors
*        of A (A = U*S*(V**H), m.ge.n)
*
         CALL ZGESVD('N','N',M,N,A,
     +               LDA,S,DUMMY,1,DUMMY,1,WORK,-1,RWORK,INFO)
      ELSE
         WRITE (NOUT,*) 'MMAX and/or NMAX too small'
      END IF
      END


and the input matrix is (file: test_zgesvd.d):
Code: Select all

   6            12                                       :Values of M and N

 (-1.377823724188512,0) (0,0) (0,0) (1,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0)
 (0,0) (-1.377823724188512,0) (0,0) (0,0) (1,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0)
 (0,0) (0,0) (-1.377823724188512,0) (0,0) (0,0) (1,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0)
 (0.01588062702035706,0) (0,0) (0,0) (-1.382229846126757,0) (0,0) (0,0) (0.005286842551427753,0) (-0.003762200585945338,0) (-0.007151269102833426,0) (0.004018458066065118,0) (-0.002859598171060272,0) (-0.005435583664416818,0)
 (0,0) (0.07630229180736924,0) (0,0) (0,0) (-1.28781829905337,0) (0,0) (0.01156945407605109,0) (-0.008233006086446525,0) (-0.01564946916158397,0) (0.04498579155453522,0) (-0.03201259914577823,0) (-0.0608502141081553,0)
 (0,0) (0,0) (0.4070646177524311,0) (0,0) (0,0) (-0.2954402733863342,0) (-0.5243179610007159,0) (0.373112934782985,0) (0.7092208247345488,0) (0.02138338930221806,0) (-0.01521675725723726,0) (-0.02892432860319349,0) :End of matrix A


The output I get is the error:
Code: Select all
 ** On entry to ZUNGLQ parameter number  5 had an illegal value

which is issued by the call to ZUNGLQ done at line 476 of ZGESVD.f (i.e., when calculating the space for ZUNGLQ)
I have not tried, but similar instructions are present in the other variant of ?GESVD.
I think the problem is related to the fact that since JOBVT='N' I pass LDVT as 1 which seems to be legal since the ?GESVD doc claims that LDVT>=1 and if JOBVT='N' the VT is not dereferenced.

A possible workaround that seems to work for me is to replace lines 476,477:
Code: Select all
            CALL ZUNGLQ( N, N, M, VT, LDVT, DUM(1), DUM(1), -1, IERR )
            LWORK_ZUNGLQ_N=DUM(1)

with:
Code: Select all
            IF( .NOT.WNTVN ) THEN
                CALL ZUNGLQ( N, N, M, VT, LDVT, DUM(1), DUM(1), -1,
     $                       IERR )
            END IF


What do you think?

Thank you very much,

-- Marco
sguazt
 
Posts: 5
Joined: Mon Aug 02, 2010 9:24 am
Location: Alessandria, Italy

Re: LAPACK trunk; issue with ?GESVD

Postby admin » Fri Nov 18, 2011 5:04 pm

Indeed, there is a problem

best would be to replace:
Code: Select all
CALL ZUNGLQ( N, N, M, VT, LDVT, DUM(1), DUM(1), -1, IERR )

by
Code: Select all
CALL ZUNGLQ( N, N, M, A, N, DUM(1), DUM(1), -1, IERR )


Thank you Marco for reporting the bug, this will be corrected very soon
Julie
admin
Site Admin
 
Posts: 504
Joined: Wed Dec 08, 2004 7:07 pm


Return to User Discussion

Who is online

Users browsing this forum: Bing [Bot] and 2 guests