LAPACK Archives

[Lapack] Bug report: LAPACK functions ZTGSEN and ZTGSYL

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




<Prev in Thread] Current Thread [Next in Thread>


For additional information you may use the LAPACK/ScaLAPACK Forum.
Or one of the mailing lists, or