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 = 1e32. 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 leastsquares solution to this matrix
1 0 1
1 1e20 0
0 0 1
and this righthand 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 = 1e32 is
chosen. I checked that the singular values of the matrix above are
1.73205 1. 5.77350e21
so this value of RCOND is small enough and it should not replace
5.77350e21 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 nontext attachment was scrubbed...
Name: bugindgelsd1.c
Type: text/xcsrc
Size: 3007 bytes
Desc: not available
Url :
http://lists.cs.utk.edu/private/lapack/attachments/20061108/799ecd56/bugindgelsd1.c
