Description
In DGEQP3, line 335-337, we do the panel factorization by JB columns among column J:N.
As noticed in DLAQPS, argument F dimension is (LDF, NB)=(N-J+1, JB), so in DGEQP3, we would expect the F starting from WORK( 2*N+JB+1 ) to WORK(LWORK) where LWORK can be the minimal value 3N+1.
But when we compute NB at line 301
NB = ( LWORK-2*SN ) / ( SN+1 )
we use SN as the free columns.
Here is a simple case, M=N=376. LWORK=3N+1 = 1129.
And input arg JPVT has the first 188 elements as non-zero, and the second 188 elements are all zero.
Then at line 170, we have NXFD=188.
So we get SM=SN=188.
Then NB = (1129-2*188)/(188+1) = 3
This will cause memory overwritten as F ranges from WORK( 2*N+JB+1 ) to WORK(LWORK) which are WORK(756) to WORK(1129). But its dimension is (LDF, NB) which is (188, 3). So we expect 564 elements but actually we can only write to 374 elements.
A quick fix is to set NB as
NB = ( LWORK-2N ) / ( SN+1 )
Here we excludes the first 2N in LWORK instead of 2SN. I verified this quick fix and resolves the problem.
A more complicate fix may be redesign the structure of workspace. So WORK(J) becomes WORK(1), WORK(N+J) becomes to WORK(SN+1) and etc. Then we only handle SN columns in the workspace.