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
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.
-------------- next part --------------
A non-text attachment was scrubbed...
Size: 3007 bytes
Desc: not available