DGEQPF and DGEQP3 from netlib's LAPACK and ATLAS

[jgpallero]

Hello,
I'm starting with LAPACK and I have some doubts about the behaviour of some functions. Attached I send a *.zip file that contains one test source code (in C) that tests DGEQRP and DGEQR3 subroutines. The problem is that if I link it with the reference LAPACK 3.2.2 and BLAS from netlib I obtain different results that if I link it against the ATLAS distribution from the repositories of my Debian GNU/Linux. The only differences in the results are about the signs of the numbers in Q and R matrices: they are different for reference LAPACK and ATLAS. The permutation vectors are the same and with both results the original matrix can be reconstructed correctly. Is the QR factorization not unique? Is my code wrong? If I try the decomposition with GNU Octave I obtain the same results that ATLAS (of course, Octave is linked against ATLAS). I compiled the reference LAPACK in my machine (iBook G4 with Debian GNU/Linux) with the parameters by default for gfortran (make.inc.gfortran), that is -O2 optimization (and -O0 when necesary).

On the other hand (for this I don't have code now), yesterday I test the function DSYEV for eigendecomposition and I obtain different eigenvectors (again only the signs was changed) but the same eigenvalues if I call DSYEV working with UPPER or LOWER part of the symmetric matrix. Is this behaviour right?

Thanks.

qrp.zip

Program for testing DGEQRP and DGEQP3 and file results

(1.94 KiB) Downloaded 49 times

[sven]

If A is real and has full column rank with m >= n, then the upper triangular matrix R of the QR factorization is unique up to sign (which is what you are observing). See theorem 5.2.2 in G H Golub and C F Van Loan, Matrix Computations, 3rd Edition, John Hopkins.

For eigenvalue problems remember that eigenvectors represent directions so again they are not unique. If

A x = lambda x then A (alpha)x = lambda (alpha)x, alpha /= 0,

so that if x is an eigenvector, then (alpha)x is also an eigenvector.

Hope that helps,

Sven Hammarling.

[jgpallero]

Sven, thank you for your reply and explanation.
I've compiled again lapack versions 3.2.1 and 3.2.2 and run my example. The absolute values are the same, but the signs are the opposite. Are the algorithms changed? The results are the same as you explained, but I don't understand the numeric behavior of the code to produce these results. Can be related with the fact that my program is a C program linking fortran code. The compiler was the same in building lapack (gfortran 4.4.4) and the C program (gcc 4.4.4).

PS: I'll upload again the files, because I can't download the originals in the first post.