Dear LAPACK team,
I would like to report bugs in the LAPACK functions ZTGSEN and ZTGSYL.
Affected are the documentation of ZTGSEN and ZTGSYL and the workspace
query in ZTGSEN.
My observations are based on the LAPACK source as available today, 20
Sep 2006, on www.netlib.org:
--------------------
SUBROUTINE ZTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
$ ALPHA, BETA, Q, LDQ, Z, LDZ, M, PL, PR, DIF,
$ WORK, LWORK, IWORK, LIWORK, INFO )
*
* -- LAPACK routine (version 3.0) --
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
* Courant Institute, Argonne National Lab, and Rice University
* June 30, 1999
--------------------
SUBROUTINE ZTGSYL( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
$ LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK,
$ IWORK, INFO )
*
* -- LAPACK routine (version 3.0) --
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
* Courant Institute, Argonne National Lab, and Rice University
* June 30, 1999
--------------------
According to http://www.netlib.org/lapack/release_notes.html the two
affected functions have not been changed since June 30, 1999.
In the following I describe the bugs:
* ZTGSYL: (http://www.netlib.org/lapack/complex16/ztgsyl.f)
---------
In the documentation to the function it says:
* IWORK (workspace) INTEGER array, dimension (M+N+2)
* If IJOB = 0, IWORK is not referenced.
however, even if IJOB = 0, ZTGSYL uses IWORK, as for example in this
code snippet (ll 325-342)
* Determine block structure of A
*
P = 0
I = 1
40 CONTINUE
IF( I.GT.M )
$ GO TO 50
P = P + 1
IWORK( P ) = I
I = I + MB
IF( I.GE.M )
$ GO TO 50
GO TO 40
50 CONTINUE
IWORK( P+1 ) = M + 1
IF( IWORK( P ).EQ.IWORK( P+1 ) )
$ P = P - 1
*
which is executed regardless off the value of IJOB. This means that
providing a too small IWORK (as, according to the documentation, it is
not referenced) can result in crahes due to overwritten memory.
* ZTGSEN: (http://www.netlib.org/lapack/complex16/ztgsen.f)
---------
Here is a similar bug to the one in ZTGSYL, but here it also affects the
workspace query. In the documentation it says:
* WORK (workspace/output) COMPLEX*16 array, dimension (LWORK)
* IF IJOB = 0, WORK is not referenced. Otherwise,
* on exit, if INFO = 0, WORK(1) returns the optimal LWORK.
and
* IWORK (workspace/output) INTEGER, dimension (LIWORK)
* IF IJOB = 0, IWORK is not referenced. Otherwise,
* on exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
Again, even if IJOB = 0, WORK and IWORK are used within the code. For
IJOB = 0, a workspace query return LWORK = 1 and LIWORK = 1 as the
optimal workspace size (ll 395 - 404):
IF( IJOB.EQ.1 .OR. IJOB.EQ.2 .OR. IJOB.EQ.4 ) THEN
LWMIN = MAX( 1, 2*M*( N-M ) )
LIWMIN = MAX( 1, N+2 )
ELSE IF( IJOB.EQ.3 .OR. IJOB.EQ.5 ) THEN
LWMIN = MAX( 1, 4*M*( N-M ) )
LIWMIN = MAX( 1, 2*M*( N-M ), N+2 )
ELSE
LWMIN = 1 <---------------
LIWMIN = 1
END IF
But even for IJOB = 0, WORK(1) and WORK(2) are referenced
as shown in this code snippet: (ll 617 - 623)
DO 60 K = 1, N
DSCALE = ABS( B( K, K ) )
IF( DSCALE.GT.SAFMIN ) THEN
WORK( 1 ) = DCONJG( B( K, K ) / DSCALE ) <----------
WORK( 2 ) = B( K, K ) / DSCALE <----------
B( K, K ) = DSCALE
CALL ZSCAL( N-K, WORK( 1 ), B( K, K+1 ), LDB )
Again, providing a Workspace of size 1 only will result in overwritten
memory and possibly a crash.
-------------------------------------------------
If you have further questions in this matter, don't hesitate to contact me.
Kind regards
Michael Wimmer
|