I currently trying to cheaply compute a good rank estimate for a matrix $A$. Therefore I compute a columnt pivoting QR decompostion using
[Q,R,E]=qr(A)
in Matlab. I estimate the rank of $A$ using
tol = size(A,n)*eps*norm(A,'fro');
r = sum(abs(diag(R))>tol)
This works fine and a plot over all diagonal entries of R looks like:
If port the whole algorithm to C/Fortran I replace [Q,R,E]=qr(A) using DGEQP3 from LAPACK, which also computes a column pivoting QR decomposition. But if I use the same estimate for the rank I mostly get something wrong. The same plot for the $R$ produced from DGEQP3 looks like
The input matrix is exactly the same for both experiments.
My question now is on which LAPACK function does the column pivoting QR decomposition from Matlab rely on?
Thanks for any help, Grisu
Edit: DGEQPF gives the same wrong result.
Edit2:
- The input matrix $A$ is dense and build as $E+sign(E,F)$
- $A$ is availble here: http://www-e.uni-magdeburg.de/makoehle/A.mtx.gz (MatrixMarket Format)
- The wrong $R$: http://www-e.uni-magdeburg.de/makoehle/R_wrong.mtx.gz
- I used LAPACK 3.4.0 with OpenBlas/GotoBLAS (64 bit)
- Matlab 7, 2007b, 2010b Linux 32bit
Edit3: - Using GDB I found out, that Matlab 2010b calls DGEQP3: #3 0xaa46ce2f in dgeqp3_ () from /usr/ubuntu10.04/matlabr2010b/bin/glnx86/../../bin/glnx86/../../bin/glnx86/mllapack.so Why do I get the wrong result using LAPACK(3.4.0 inlcude the fixes mentioned in Working Note 176)?