LAPACK/ScaLAPACK Development

Hi,

I had already reported this kind of bug (for DGELSD) by mail a few months ago, but it could not be reproduced (by Julien Langou), and I did not have time to follow up.

Since I have run across it again with DGESDD, I have investigated it more closely, and am still able to reproduce it with DGELSD, too:

When an optimum workspace query is requested by passing "-1" as argument for LWORK, both functions return sizes that are often smaller than the minimum workspace size as defined by the formulas in the documentation. I am able to reproduce this on Linux (i686) with the latest LAPACK-rpm for Fedora Core 3 (lapack-3.0-28).

Here is example code for DGELSD:

Code: Select all

      double precision WORK(1)
      LDA = 100
      LDB = 100
      M = 100
      N = 1
      NRHS = 1
      LWORK = -1
      CALL DGELSD(M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK,
     &            LWORK, IWORK, INFO)
      WRITE(*,*) 'M=',M,', N=',N,', NRHS=',NRHS
      WRITE(*,*) 'IF LWORK = -1 THEN WORK(1) = ', WORK(1)
      end


The above will print 67 as optimum workspace size. But the minimum length according to the documentation should be 739.

Here is example code for DGESDD:

Code: Select all

      double precision WORK(1)
      LDA = 50
      M = 50
      N = 100
      LDU = 50
      LDVT = 100
      LWORK = -1
      CALL DGESDD('A', M, N, A, LDA, S, U, LDU, VT, LDVT, WORK,
     &            LWORK, IWORK, INFO)
      WRITE(*,*) 'M=',M,', N=',N
      WRITE(*,*) 'IF LWORK = -1 THEN WORK(1) = ', WORK(1)
      end


Here it prints 10450 as optimum workspace size. But the formula in the documentation for the minimum demands 17700.

Has anybody else seen these problems? Could one of the developers please investigate it? - Thanks!

Best regards,
Markus

Markus,

You are using the latest version of Lapack with the bug fixes from http://www.netlib.org/lapack/release_notes.htm.
As you can see on the release notes, both of the routines, DGELSD and DGELSD have been corrected regarding to a error in LQUERY and setting of WORK(1).
Here is the code used in your library: http://www.netlib.org/lapack/patch/src/dgesdd.f and http://www.netlib.org/lapack/patch/src/dgelsd.f

Unfortunately, the documentation (Lapack User Guide) has not been updated, neither the man pages. The only way to know the formula of the workspace is to look directly at the source file.

Here are the results I obtain with your examples

Code: Select all

 TEST DGELSD
 M= 100, N= 1, NRHS= 1
 IF LWORK = -1 THEN WORK(1) =   67.

  ************ OLD FORMULA *************
 12*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + (SMLSIZ+1)**2
 ******* OLD OPT_SPACE =  739 *******
 
 ************ NEW FORMULA *************
 If M >= N, LWORK >= 11*N+2*N*SMLSIZ+8*N*NLVL+N*NRHS
 ******* NEW OPT_SPACE =  62 *******



Code: Select all

 TEST DGESDD
 M= 50, N= 100
 IF LWORK = -1 THEN WORK(1) =   10450.

  ************ OLD FORMULA *************
 JOBZ = A so
 LWORK >= 3*min(M,N)*min(M,N) +
       max(max(M,N),4*min(M,N)*min(M,N)+4*min(M,N)).
 ******* OLD OPT_SPACE =  17550 *******

  ************ NEW FORMULA *************
 LWORK >= 4*min(M,N)*min(M,N) + max(M,N) + 9*min(M,N)
 ******* NEW OPT_SPACE =  10550 *******


We are currently working on a future release of Lapack, everything will be updated.

Sorry for the inconvenience
Juliie