Out of bounds array access in ZGESVD for 1x1 matrix

[frankotto]

Hello,

if LAPACK is compiled with array boundary checks (e.g. gfortran option -fcheck=bounds), then ZGESVD will cause an out-of-bounds array access when applied to a 1x1 matrix. (I'm requesting the right but not the left singular vectors.) Though this is an admittedly trivial (and maybe pointless) operation, nothing in the documentation points out that the matrix must be larger than 1x1. And in fact, in the interest of generality, the routine should still do the right thing even for the 1x1 case.

Example code follows (I would like to attach it but .f90 is not allowed...)

program svdtest
complex*16 :: A(1,1), U(1), V(1), LWORKC(1)
real*8 :: S(1), RWORK(5)
integer :: M, N, LDA, LDU, LDV, LWORK, INFO
complex*16,allocatable :: WORK(:)
M = 1
N = 1
LDA = 1
LDU = 1
LDV = 1
A(1,1) = (42.d0, 0.d0)
print *, "A =", A
print *, "workspace query..."
call ZGESVD('N', 'O', M, N, A, LDA, S, U, LDU, V, LDV, LWORKC, -1, RWORK, INFO)
print *, "INFO =", INFO
LWORK = max(1,int(LWORKC(1)))
print *, "LWORK =", LWORK
allocate(WORK(LWORK))
print *, "ZGESVD..."
call ZGESVD('N', 'O', M, N, A, LDA, S, U, LDU, V, LDV, WORK, LWORK, RWORK, INFO)
print *, "INFO =", INFO
print *, "S =", S
print *, "A =", A
deallocate(WORK)
end program svdtest

LAPACK is compiled with "OPTS = -O1 -frecursive -fcheck=bounds" in make.inc, and uses the refblas.
Tested with version 3.6.0 of LAPACK.
gfortran is version 4.8.4.

Test program compilation:
gfortran -Wall -fcheck=bounds svdtest.f90 lapack-3.6.0/liblapack.a lapack-3.6.0/librefblas.a

Test program invocation:
A = ( 42.0000000000000000 , 0.0000000000000000 )
workspace query...
INFO = 0
LWORK = 66
ZGESVD...
At line 683 of file zgesvd.f
Fortran runtime error: Index '2' of dimension 1 of array 'a' above upper bound of 1

And indeed, line 683 of zgesvd.f reads:
CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ), LDA )

Expected behavior:
LAPACK should not reference array entries other than A(1,1).

Thank you and best regards,
Frank

[lawrence mulholland]

Yes, there are several places in each of [cdsz]gesvd.f
where calls with args n-1,n-1,.. a(2,1) are made.

some are already protected with If (n .gt. 1)

the remainder need to be also.

I have made these changes and have committed them as revision r1685

Lawrence (NAG)

[Julien Langou]

Thanks Lawrence. This is great! Julien.