LAPACK Archives

[Lapack] A bug in dgelsd/dlalsd?

Dear LAPACK authors,

I think that there is an error in the LAPACK function dgelsd (to be more
precise, in dlalsd called from dgelsd).

The first problem I encountered is that the argument RCOND cannot be
zero. Here is the dgelsd prototype:

SUBROUTINE DGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
      $                   WORK, LWORK, IWORK, INFO )

In the documentation of dgelsd it is written that RCOND can be zero.

*  RCOND   (input) DOUBLE PRECISION
*          RCOND is used to determine the effective rank of A.
*          Singular values S(i) <= RCOND*S(1) are treated as zero.
*          If RCOND < 0, machine precision is used instead.

But in the code if RCOND is zero it is immediately replaced by EPS which 
is not documented. This replacement is done in dlalsd:

*
*     Set up the tolerance.
*
       IF( ( RCOND.LE.ZERO ) .OR. ( RCOND.GE.ONE ) ) THEN
          RCOND = EPS
       END IF

I think that RCOND.LT.ZERO should be used above.

But since I cannot change this code (I use the MKL library from Intel in
binary form) and I find some cases where RCOND = 0 is useful, I use a 
small number for such problems, say, RCOND = 1e-32. But even with this 
small value dgelsd still does not work as I expected. I am sending you a 
standalone C program I wrote which shows the problem. (I usually write 
programs in C, not FORTRAN, so it was easier for me to write it in C).

In this program I want to find the least-squares solution to this matrix

1  0      1
1  1e-20  0
0  0      1

and this right-hand side

0 1 0

This is a simplified example of a much larger system I wanted to
solve. Solving this problem by hand it is easy to find the exact solution

0 1e+20 0

So I expected a similar solution from dgelsd when RCOND = 1e-32 is 
chosen. I checked that the singular values of the matrix above are

  1.73205  1.  5.77350e-21

so this value of RCOND is small enough and it should not replace 
5.77350e-21 with zero. But the solution given by dgelsd is always the same

0.666667  0  -0.333333

independently of what value of RCOND I choose. It also replaces the 
smallest singular value by zero. And to be sure that it is not a bug in 
MKL I linked the standalone program with the original CLAPACK (ver. 3.0) 
library. The result is still the same.

I suspect that the problem is in the function dlalsd where in many 
places EPS (= DLAMCH( 'Epsilon' )) is used instead of TOL (but I can be 
wrong here). TOL is defined in dlalsd as

TOL = RCOND*ABS( D( IDAMAX( N, D, 1 ) ) )

If it is a known problem, please let me know. If it is not known, I can 
provide more details of why I think it is a buggy behavior and I can 
look at the code in dlalsd to show which lines I think are not correct.

Best regards,
Zbigniew Leyk
-------------- next part --------------
A non-text attachment was scrubbed...
Name: bug-in-dgelsd1.c
Type: text/x-csrc
Size: 3007 bytes
Desc: not available
Url : 
http://lists.cs.utk.edu/private/lapack/attachments/20061108/799ecd56/bug-in-dgelsd1.c
 

<Prev in Thread] Current Thread [Next in Thread>
  • [Lapack] A bug in dgelsd/dlalsd?, Zbigniew Leyk <=


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