Bug (maybe) in xGESVD doc: wrong minimum workspace size

[sguazt]

Dear all,

It seems that the minimum workspace size in xGESVD doc is wrong.
E.g., in DGESVD I read:

Code: Select all

00116 *  LWORK   (input) INTEGER
00117 *          The dimension of the array WORK.
00118 *          LWORK >= MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N)).


So for an input matrix with M=1 and N=3, LWORK should be >= 6.
However, DGESVD, if called with LWORK=-1, tell me that the optimal size is 5 (WORK(1) == 5).

Here below there is a sample code test_dgesvd.f ...

Code: Select all

      PROGRAM TEST_DGESVD
*     DGESVD Example Program Text
*     .. Parameters ..
      CHARACTER        JOBU, JOBVT
      PARAMETER        (JOBU='N',JOBVT='N')
      INTEGER          NIN, NOUT
      PARAMETER        (NIN=5,NOUT=6)
      INTEGER          MMAX, NMAX
      PARAMETER        (MMAX=10,NMAX=8)
      INTEGER          LDA, LDU, LDVT
      PARAMETER        (LDA=MMAX,LDU=MMAX,LDVT=NMAX)
*     .. Local Scalars ..
      INTEGER          I, INFO, J, LWKOPT, LWKMIN, M, N
*     .. Local Arrays ..
      DOUBLE PRECISION A(LDA,NMAX), U(LDU,MMAX), S(NMAX),
     +                 VT(LDVT,NMAX), WORK(MMAX*NMAX)
*     .. External Subroutines ..
      EXTERNAL         DGESVD
*     .. Executable Statements ..
      WRITE (NOUT,*) 'DGESVD Example Program Results'
      WRITE (NOUT,*)
*     Skip heading in data file
      READ (NIN,*)
      READ (NIN,*) M, 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**T))
*
         LWKOPT=-1
         CALL DGESVD(JOBU,JOBVT,M,N,A,
     +               LDA,S,U,LDU,VT,LDVT,WORK,LWKOPT,INFO)
         LWKMIN= MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N))
         LWKOPT=WORK(1)
         CALL DGESVD(JOBU,JOBVT,M,N,A,
     +               LDA,S,U,LDU,VT,LDVT,WORK,LWKOPT,INFO)
*
         IF (INFO.EQ.0) THEN
*
*           Print solution
*
            WRITE (NOUT,*) 'Singular values'
            WRITE (NOUT,99999) (S(J),J=1,N)
*
            IF (JOBU.EQ.'A') THEN
                WRITE (NOUT,*) 'Left Singular vectors'
                DO I=1,M
                    WRITE (NOUT,99999) (U(I,J),J=1,M)
                END DO
            ELSE IF (JOBU.EQ.'S') THEN
                WRITE (NOUT,*) 'Left Singular vectors'
                DO I=1,M
                    WRITE (NOUT,99999) (U(I,J),J=1,MIN(M,N))
                END DO
            END IF

            IF (JOBVT.EQ.'A') THEN
                WRITE (NOUT,*) 'Right Singular vectors'
                DO I=1,N
                    WRITE (NOUT,99999) (VT(I,J),J=1,N)
                END DO
            ELSE IF (JOBVT.EQ.'S') THEN
                WRITE (NOUT,*) 'Right Singular vectors'
                DO I=1,MIN(M,N)
                    WRITE (NOUT,99999) (VT(I,J),J=1,N)
                END DO
            END IF
         ELSE
            WRITE (NOUT,99998) 'Failure in DGESVD. INFO =', INFO
         END IF
*
*        Print workspace information
*
         IF (LWKOPT.LT.LWKMIN) THEN
            WRITE (NOUT,*)
            WRITE (NOUT,99997) 'Optimum workspace required = ', LWKOPT,
     +                         'Minimum Workspace expected = ', LWKMIN
         END IF
      ELSE
         WRITE (NOUT,*) 'MMAX and/or NMAX too small'
      END IF
      STOP
*
99999 FORMAT (3X,(8F8.4))
99998 FORMAT (1X,A,I4)
99997 FORMAT (1X,A,I5,/1X,A,I5)
      END


... and an input file test_dgesvd.data

Code: Select all

DGESVD Example Program Data

   1      3            :Values of M and N

  -1.00   0.00   1.00  :End of matrix A


Run as:

Code: Select all

$ ./test_dgesvd < test_dgesvd.data


(tested with gfortran 4,4,4 and LAPACK 3.2.1)

What do you think?

Thank you very much!!

Best,

-- Marco

[admin]

Hi Marco,
I couldn't reproduce your problem (I believe you changed the value of NB in ILAENV) but I saw where the mistake comes from.
The formula LWORK >= MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N)). in the comment is very "conservative".
It merged all the cases. For example for your case (Path 1t :* Path 1t(N much larger than M, JOBVT='N') ), the minimum workspace is given by

Code: Select all

MINWRK = MAX( 4*M, BDSPAC ) where BDSPAC = 5*M

this is where the 5 is coming from.
In general you have to use the LWORK query like you did. Comments can be incomplete.
I will add a little note in the comments saying that this formula is a conservative minimum worksize.
This will be listed as bug0077 in our LAPACK Errata.
Sorry for answering so late but I didn't understand your problem in the first place.
Julie

[sguazt]

Thank you Julie!

Sorry for the ill-posed question.
The reason for this question is that actually I am using LAPACK via a third-party C++ library which performs some check on the minimum workspace size.
Such checks seem to relies on the size requirements found in the comments.

I am trying to contact the author of the library in order to fix this.

Thank you very much for the feedback.

Best,

-- Marco