DGELSD & DGESDD - bug in optimum workspace computation?

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

DGELSD & DGESDD - bug in optimum workspace computation?

Postby markusm73 » Wed Jun 15, 2005 11:35 am

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
markusm73
 
Posts: 8
Joined: Tue Jun 14, 2005 8:10 pm
Location: New York

Postby Julie » Mon Jun 20, 2005 11:18 am

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
Julie
 
Posts: 299
Joined: Wed Feb 23, 2005 12:32 am
Location: ICL, Denver. Colorado


Return to User Discussion

Who is online

Users browsing this forum: Google [Bot] and 1 guest

cron