The problem I encountered is that I cannot pass the argument RCOND as zero in DGELSD. But in the documentation of DGELSD it is written that RCOND can be zero.

- Code: Select all
`* 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.

Here is the DGELSD prototype:

- Code: Select all
`SUBROUTINE DGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,`

$ WORK, LWORK, IWORK, INFO )

In the DGELSD code, if RCOND is zero it is immediately replaced by EPS and this behavior is not documented. This replacement is done in DLALSD.

- Code: Select all
`*`

* Set up the tolerance.

*

IF( ( RCOND.LE.ZERO ) .OR. ( RCOND.GE.ONE ) ) THEN

RCOND = EPS

END IF

I think that in place of

- Code: Select all
`RCOND.LE.ZERO`

- Code: Select all
`RCOND.LT.ZERO`

This small bug implies another problem. Since I cannot change LAPACK code (I use the MKL library from Intel in binary form) and I find some cases where RCOND = 0 is useful for me, I try to assign a small number to RCOND, say, RCOND = 1e-32. But even with this small value DGELSD still does not work as expected.

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

- Code: Select all

[ 1 0 1]

A = [ 1 1e-20 0]

[ 0 0 1]

and to this right-hand side

- Code: Select all
`b = [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 that the exact solution is

- Code: Select all
`x = [0 1e+20 0]`

So I expected a similar solution from DGELSD with RCOND = 1e-32. I checked that the singular values of the matrix A are

- Code: Select all
`s = [1.73205 1. 5.77350e-21]`

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

- Code: Select all
`x1 = [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 my 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

- Code: Select all
`TOL = RCOND*ABS( D( IDAMAX( N, D, 1 ) ) )`

I am interested if this is a known problem.

Zbigniew