I'm trying to compute the solution of a singular square system using a QR factorization. When I'm using a QR factorization with pivoting (i.e. ?geqp3), everything works OK. However, I'm concerned about the cost of pivoting, and I'd rather use ?geqrf. The problem is that the computed vector is in this case not a solution of the original system, as shown in the code attached. Do you see what I'm doing wrong ?
Also, could you comment on the performance of ?geqp3 vs. ?geqrf ?
The code in the example operates with the matrix like it is stored in row-major storage, however LAPACK interfaces assume matrices are in column-major storage. (Some details could found at https://software.intel.com/en-us/node/520868, Matrix Layout paragraph )
Replacing row-major "A = A = 0;" with column-major "A = A = 0;" (and adding "x = 1.0;") results correct solution.
You may want to try LAPACKE interface if you would like to operate with input matrices stored row-major.
As of QRF and QP3: QRF is faster. Difference could be up to 2x in some cornercases (too small matrix, or medium matrix and many threads). Introducing fixed (JPVT(J).ne.0) columns in QP3 factorization improves it performance: the more fixed columns, the closer QP3 performance to QR.