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 ```
 Current 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