potential memory over written in GEQP3

[winice-test]

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.

[winice-test]

A simple case to show overwritten. We allocate slightly larger workspace as 3N+1+8. Then we found the tail 8 elements are changed after calling dgeqp3.

  program xdgeqp3
  implicit  none
  real*8 A(320,320)
  real*8 TAU(320)
  integer i, n, lwork, info
  integer jpvt(320)
  real*8 work(961+8)

  n = 320
  work(961+1:961+8) = 1.D0
  lwork = 961
  JPVT(1:n) = 0
  JPVT(1:n/2) = 1

  A = 0.0D0
  DO i=1, n
    A(i,i) = 1.0D0
  ENDDO

  write(6,'(A20,8F6.2)')'orig tail = ', work(961+1:961+8)
  call dgeqp3 (n,n,A,n,JPVT,TAU,work,lwork,INFO)
  write(6,*)'INFO = ', INFO
  write(6,'(A20,8F6.2)')'new  tail = ', work(961+1:961+8)

  end program

Output on Ubuntu1604 with blas-3.6.0-2ubuntu2, lapack-3.6.0-2ubuntu2 and gfortran-5-5.4.0-6ubuntu1~16.04.11.

guest@localhost:~/tmp$ gfortran -o xqp3 xqp3.f -llapack -lblas guest@localhost:~/tmp$ ./xqp3 orig tail = 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 INFO = 0 new tail = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00