LAPACK/ScaLAPACK Development

I need to solve some linear least squares problems. The matrix is not necessarily full rank, so I guess I need to use xGELSY or xGELSD.

For xGELSY, I don't know what to pass for the JPVT and RCOND input parameters.

For xGELSD, I don't know what to pass for RCOND or what dimension to use for IWORK. How come the function can compute the size of WORK for you, but not IWORK?

I see that
LIWORK >= 3 * MINMN * NLVL + 11 * MINMN
where
MINMN = MIN( M,N)
and
NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) +1)
but what is SMLSIZ?

By looking at the DGELSY example at http://www.nag.com/lapack-ex/node48.html, I can see that I can initialize JPVT to all zeros. I'm still not sure about RCOND. The example uses RCOND = 0.01 and comes up with the solution:

0.6344 0.9699 -1.4402 3.3678 3.3992

If I use RCOND = 0.01, or 0.001, I get about that answer. But if I use RCOND = 0.0001, I get

-0.799781 -3.288069 -7.475135 4.939312 0.767767

Now, I wouldn't expect to always get exactly the same solution, but since this algorithm is supposed to be finding a minimum-norm solution, I'd expect the norms to be fairly close... and they're not.

RCOND is used to determine the numerical rank of the matrix A, and should be chosen to reflect the accuracy of your data.

Consider the example where

A = ( 1 0 ), b = ( 1 ).
( 0 eps ) ( 1 )

If eps is small, but treated as nonzero, so the rank of A is taken to be 2, then the solution is

x = 1, y = 1/eps.

On the other hand, if we treat eps as zero, so that we regard the rank of A as 1, then the solution is

x = 1, y = 0.

So you see that we can get quite a difference in solution.

When I run your example, I found that the value returned in RANK changed from 4 to 5 when RCOND = 0.0001.

I hope that helps,

Sven Hammarling.

Thanks, that helps.