Correct bug reported by Alex Zotkevich from Intel

See: http://icl.utk.edu/lapack-forum/viewtopic.php?f=13&t=4392

Bug: During workspace computation, LAPACK code CGESVD was calling other
subroutines (e.g. CGEQRF) with REAL DUM variable as COMPLEX WORK variable.  DUM
(in CGESVD) is REAL while WORK (in called subroutines) is COMPLEX.  This
corrupts the stack when a value is set in WORK.

Fix: In CGESVD, use the COMPLEX CDUM variable (already present in the code)
instead of the REAL DUM variable.  Since I was at it, the COMPLEX "TAU"
variables (not referenced anyway) were passed the REAL DUM variable, I changed
the code so that the COMPLEX CDUM variable is passed. This is cleaner like
this.

Same problem with ZGESVD. Same fix.

Alex's post:
Hi, We recently found a stack corruption issue in (C,Z)GESVD that potentially
could even lead to incorrect xerbla error message.  In ZGESVD array DUM which
is used in LWORK query is a double precision array of size 1 allocated on
stack:
	DOUBLE PRECISION DUM( 1 )
DUM comes to ( ZGEQRF, ZUNGQR, ... ) as a WORK array to return an optimal LWORK
value.  But in ( ZGEQRF, ZUNGQR, ... ) array WORK is declared as a COMPLEX*16
array. So WORK(1) = 1 corrupts the stack as it deals with complex value while
pointer on input of the function is a pointer to double: (oooooooo|xxxxxxxx),
oooooooo fills with LWORK value, xxxxxxxx corrupts. Let compiler use xxxxxxxx
to hold some value. After LWORK query the value will turn to be a zero.  "Hacky
fix" would be to allocate DUM array of size 2.
W.B.R.
Alex Zotkevich

2 files changed, 98 insertions, 96 deletions

diff --git a/SRC/cgesvd.f b/SRC/cgesvd.f
index 25b178dd..47062eb5 100644
--- a/SRC/cgesvd.f
+++ b/SRC/cgesvd.f
@@ -321,24 +321,24 @@
 *
             MNTHR = ILAENV( 6, 'CGESVD', JOBU // JOBVT, M, N, 0, 0 )
 *           Compute space needed for CGEQRF
-            CALL CGEQRF( M, N, A, LDA, DUM(1), DUM(1), -1, IERR )
-            LWORK_CGEQRF=DUM(1)
+            CALL CGEQRF( M, N, A, LDA, CDUM(1), CDUM(1), -1, IERR )
+            LWORK_CGEQRF=CDUM(1)
 *           Compute space needed for CUNGQR
-            CALL CUNGQR( M, N, N, A, LDA, DUM(1), DUM(1), -1, IERR )
-            LWORK_CUNGQR_N=DUM(1)
-            CALL CUNGQR( M, M, N, A, LDA, DUM(1), DUM(1), -1, IERR )
-            LWORK_CUNGQR_M=DUM(1)
+            CALL CUNGQR( M, N, N, A, LDA, CDUM(1), CDUM(1), -1, IERR )
+            LWORK_CUNGQR_N=CDUM(1)
+            CALL CUNGQR( M, M, N, A, LDA, CDUM(1), CDUM(1), -1, IERR )
+            LWORK_CUNGQR_M=CDUM(1)
 *           Compute space needed for CGEBRD
-            CALL CGEBRD( N, N, A, LDA, S, DUM(1), DUM(1),
-     $                   DUM(1), DUM(1), -1, IERR )
-            LWORK_CGEBRD=DUM(1)
+            CALL CGEBRD( N, N, A, LDA, S, DUM(1), CDUM(1),
+     $                   CDUM(1), CDUM(1), -1, IERR )
+            LWORK_CGEBRD=CDUM(1)
 *           Compute space needed for CUNGBR
-            CALL CUNGBR( 'P', N, N, N, A, LDA, DUM(1),
-     $                   DUM(1), -1, IERR )
-            LWORK_CUNGBR_P=DUM(1)
-            CALL CUNGBR( 'Q', N, N, N, A, LDA, DUM(1),
-     $                   DUM(1), -1, IERR )
-            LWORK_CUNGBR_Q=DUM(1)
+            CALL CUNGBR( 'P', N, N, N, A, LDA, CDUM(1),
+     $                   CDUM(1), -1, IERR )
+            LWORK_CUNGBR_P=CDUM(1)
+            CALL CUNGBR( 'Q', N, N, N, A, LDA, CDUM(1),
+     $                   CDUM(1), -1, IERR )
+            LWORK_CUNGBR_Q=CDUM(1)
 *
             MNTHR = ILAENV( 6, 'CGESVD', JOBU // JOBVT, M, N, 0, 0 )
             IF( M.GE.MNTHR ) THEN
@@ -444,20 +444,20 @@
 *
 *              Path 10 (M at least N, but not much larger)
 *
-               CALL CGEBRD( M, N, A, LDA, S, DUM(1), DUM(1),
-     $                   DUM(1), DUM(1), -1, IERR )
-               LWORK_CGEBRD=DUM(1)
+               CALL CGEBRD( M, N, A, LDA, S, DUM(1), CDUM(1),
+     $                   CDUM(1), CDUM(1), -1, IERR )
+               LWORK_CGEBRD=CDUM(1)
                MAXWRK = 2*N + LWORK_CGEBRD
                IF( WNTUS .OR. WNTUO ) THEN
-                  CALL CUNGBR( 'Q', M, N, N, A, LDA, DUM(1),
-     $                   DUM(1), -1, IERR )
-                  LWORK_CUNGBR_Q=DUM(1)
+                  CALL CUNGBR( 'Q', M, N, N, A, LDA, CDUM(1),
+     $                   CDUM(1), -1, IERR )
+                  LWORK_CUNGBR_Q=CDUM(1)
                   MAXWRK = MAX( MAXWRK, 2*N+LWORK_CUNGBR_Q )
                END IF
                IF( WNTUA ) THEN
-                  CALL CUNGBR( 'Q', M, M, N, A, LDA, DUM(1),
-     $                   DUM(1), -1, IERR )
-                  LWORK_CUNGBR_Q=DUM(1)
+                  CALL CUNGBR( 'Q', M, M, N, A, LDA, CDUM(1),
+     $                   CDUM(1), -1, IERR )
+                  LWORK_CUNGBR_Q=CDUM(1)
                   MAXWRK = MAX( MAXWRK, 2*N+LWORK_CUNGBR_Q )
                END IF
                IF( .NOT.WNTVN ) THEN
@@ -471,25 +471,26 @@
 *
             MNTHR = ILAENV( 6, 'CGESVD', JOBU // JOBVT, M, N, 0, 0 )
 *           Compute space needed for CGELQF
-            CALL CGELQF( M, N, A, LDA, DUM(1), DUM(1), -1, IERR )
-            LWORK_CGELQF=DUM(1)
+            CALL CGELQF( M, N, A, LDA, CDUM(1), CDUM(1), -1, IERR )
+            LWORK_CGELQF=CDUM(1)
 *           Compute space needed for CUNGLQ
-            CALL CUNGLQ( N, N, M, DUM(1), N, DUM(1), DUM(1), -1, IERR )
-            LWORK_CUNGLQ_N=DUM(1)
-            CALL CUNGLQ( M, N, M, A, LDA, DUM(1), DUM(1), -1, IERR )
-            LWORK_CUNGLQ_M=DUM(1)
+            CALL CUNGLQ( N, N, M, CDUM(1), N, CDUM(1), CDUM(1), -1,
+     $                   IERR )
+            LWORK_CUNGLQ_N=CDUM(1)
+            CALL CUNGLQ( M, N, M, A, LDA, CDUM(1), CDUM(1), -1, IERR )
+            LWORK_CUNGLQ_M=CDUM(1)
 *           Compute space needed for CGEBRD
-            CALL CGEBRD( M, M, A, LDA, S, DUM(1), DUM(1),
-     $                   DUM(1), DUM(1), -1, IERR )
-            LWORK_CGEBRD=DUM(1)
+            CALL CGEBRD( M, M, A, LDA, S, DUM(1), CDUM(1),
+     $                   CDUM(1), CDUM(1), -1, IERR )
+            LWORK_CGEBRD=CDUM(1)
 *            Compute space needed for CUNGBR P
-            CALL CUNGBR( 'P', M, M, M, A, N, DUM(1),
-     $                   DUM(1), -1, IERR )
-            LWORK_CUNGBR_P=DUM(1)
+            CALL CUNGBR( 'P', M, M, M, A, N, CDUM(1),
+     $                   CDUM(1), -1, IERR )
+            LWORK_CUNGBR_P=CDUM(1)
 *           Compute space needed for CUNGBR Q
-            CALL CUNGBR( 'Q', M, M, M, A, N, DUM(1),
-     $                   DUM(1), -1, IERR )
-            LWORK_CUNGBR_Q=DUM(1)
+            CALL CUNGBR( 'Q', M, M, M, A, N, CDUM(1),
+     $                   CDUM(1), -1, IERR )
+            LWORK_CUNGBR_Q=CDUM(1)
             IF( N.GE.MNTHR ) THEN
                IF( WNTVN ) THEN
 *
@@ -593,21 +594,21 @@
 *
 *              Path 10t(N greater than M, but not much larger)
 *
-               CALL CGEBRD( M, N, A, LDA, S, DUM(1), DUM(1),
-     $                   DUM(1), DUM(1), -1, IERR )
-               LWORK_CGEBRD=DUM(1)
+               CALL CGEBRD( M, N, A, LDA, S, DUM(1), CDUM(1),
+     $                   CDUM(1), CDUM(1), -1, IERR )
+               LWORK_CGEBRD=CDUM(1)
                MAXWRK = 2*M + LWORK_CGEBRD
                IF( WNTVS .OR. WNTVO ) THEN
 *                Compute space needed for CUNGBR P
-                 CALL CUNGBR( 'P', M, N, M, A, N, DUM(1),
-     $                   DUM(1), -1, IERR )
-                 LWORK_CUNGBR_P=DUM(1)
+                 CALL CUNGBR( 'P', M, N, M, A, N, CDUM(1),
+     $                   CDUM(1), -1, IERR )
+                 LWORK_CUNGBR_P=CDUM(1)
                  MAXWRK = MAX( MAXWRK, 2*M+LWORK_CUNGBR_P )
                END IF
                IF( WNTVA ) THEN
-                 CALL CUNGBR( 'P', N,  N, M, A, N, DUM(1),
-     $                   DUM(1), -1, IERR )
-                 LWORK_CUNGBR_P=DUM(1)
+                 CALL CUNGBR( 'P', N,  N, M, A, N, CDUM(1),
+     $                   CDUM(1), -1, IERR )
+                 LWORK_CUNGBR_P=CDUM(1)
                  MAXWRK = MAX( MAXWRK, 2*M+LWORK_CUNGBR_P )
                END IF
                IF( .NOT.WNTUN ) THEN
diff --git a/SRC/zgesvd.f b/SRC/zgesvd.f
index b3fd8e37..6877f559 100644
--- a/SRC/zgesvd.f
+++ b/SRC/zgesvd.f
@@ -321,24 +321,24 @@
 *
             MNTHR = ILAENV( 6, 'ZGESVD', JOBU // JOBVT, M, N, 0, 0 )
 *           Compute space needed for ZGEQRF
-            CALL ZGEQRF( M, N, A, LDA, DUM(1), DUM(1), -1, IERR )
-            LWORK_ZGEQRF=DUM(1)
+            CALL ZGEQRF( M, N, A, LDA, CDUM(1), CDUM(1), -1, IERR )
+            LWORK_ZGEQRF=CDUM(1)
 *           Compute space needed for ZUNGQR
-            CALL ZUNGQR( M, N, N, A, LDA, DUM(1), DUM(1), -1, IERR )
-            LWORK_ZUNGQR_N=DUM(1)
-            CALL ZUNGQR( M, M, N, A, LDA, DUM(1), DUM(1), -1, IERR )
-            LWORK_ZUNGQR_M=DUM(1)
+            CALL ZUNGQR( M, N, N, A, LDA, CDUM(1), CDUM(1), -1, IERR )
+            LWORK_ZUNGQR_N=CDUM(1)
+            CALL ZUNGQR( M, M, N, A, LDA, CDUM(1), CDUM(1), -1, IERR )
+            LWORK_ZUNGQR_M=CDUM(1)
 *           Compute space needed for ZGEBRD
-            CALL ZGEBRD( N, N, A, LDA, S, DUM(1), DUM(1),
-     $                   DUM(1), DUM(1), -1, IERR )
-            LWORK_ZGEBRD=DUM(1)
+            CALL ZGEBRD( N, N, A, LDA, S, DUM(1), CDUM(1),
+     $                   CDUM(1), CDUM(1), -1, IERR )
+            LWORK_ZGEBRD=CDUM(1)
 *           Compute space needed for ZUNGBR
-            CALL ZUNGBR( 'P', N, N, N, A, LDA, DUM(1),
-     $                   DUM(1), -1, IERR )
-            LWORK_ZUNGBR_P=DUM(1)
-            CALL ZUNGBR( 'Q', N, N, N, A, LDA, DUM(1),
-     $                   DUM(1), -1, IERR )
-            LWORK_ZUNGBR_Q=DUM(1)
+            CALL ZUNGBR( 'P', N, N, N, A, LDA, CDUM(1),
+     $                   CDUM(1), -1, IERR )
+            LWORK_ZUNGBR_P=CDUM(1)
+            CALL ZUNGBR( 'Q', N, N, N, A, LDA, CDUM(1),
+     $                   CDUM(1), -1, IERR )
+            LWORK_ZUNGBR_Q=CDUM(1)
 *
             IF( M.GE.MNTHR ) THEN
                IF( WNTUN ) THEN
@@ -443,20 +443,20 @@
 *
 *              Path 10 (M at least N, but not much larger)
 *
-               CALL ZGEBRD( M, N, A, LDA, S, DUM(1), DUM(1),
-     $                   DUM(1), DUM(1), -1, IERR )
-               LWORK_ZGEBRD=DUM(1)
+               CALL ZGEBRD( M, N, A, LDA, S, DUM(1), CDUM(1),
+     $                   CDUM(1), CDUM(1), -1, IERR )
+               LWORK_ZGEBRD=CDUM(1)
                MAXWRK = 2*N + LWORK_ZGEBRD
                IF( WNTUS .OR. WNTUO ) THEN
-                  CALL ZUNGBR( 'Q', M, N, N, A, LDA, DUM(1),
-     $                   DUM(1), -1, IERR )
-                  LWORK_ZUNGBR_Q=DUM(1)
+                  CALL ZUNGBR( 'Q', M, N, N, A, LDA, CDUM(1),
+     $                   CDUM(1), -1, IERR )
+                  LWORK_ZUNGBR_Q=CDUM(1)
                   MAXWRK = MAX( MAXWRK, 2*N+LWORK_ZUNGBR_Q )
                END IF
                IF( WNTUA ) THEN
-                  CALL ZUNGBR( 'Q', M, M, N, A, LDA, DUM(1),
-     $                   DUM(1), -1, IERR )
-                  LWORK_ZUNGBR_Q=DUM(1)
+                  CALL ZUNGBR( 'Q', M, M, N, A, LDA, CDUM(1),
+     $                   CDUM(1), -1, IERR )
+                  LWORK_ZUNGBR_Q=CDUM(1)
                   MAXWRK = MAX( MAXWRK, 2*N+LWORK_ZUNGBR_Q )
                END IF
                IF( .NOT.WNTVN ) THEN
@@ -470,25 +470,26 @@
 *
             MNTHR = ILAENV( 6, 'ZGESVD', JOBU // JOBVT, M, N, 0, 0 )
 *           Compute space needed for ZGELQF
-            CALL ZGELQF( M, N, A, LDA, DUM(1), DUM(1), -1, IERR )
-            LWORK_ZGELQF=DUM(1)
+            CALL ZGELQF( M, N, A, LDA, CDUM(1), CDUM(1), -1, IERR )
+            LWORK_ZGELQF=CDUM(1)
 *           Compute space needed for ZUNGLQ
-            CALL ZUNGLQ( N, N, M, DUM(1), N, DUM(1), DUM(1), -1, IERR )
-            LWORK_ZUNGLQ_N=DUM(1)
-            CALL ZUNGLQ( M, N, M, A, LDA, DUM(1), DUM(1), -1, IERR )
-            LWORK_ZUNGLQ_M=DUM(1)
+            CALL ZUNGLQ( N, N, M, CDUM(1), N, CDUM(1), CDUM(1), -1,
+     $                   IERR )
+            LWORK_ZUNGLQ_N=CDUM(1)
+            CALL ZUNGLQ( M, N, M, A, LDA, CDUM(1), CDUM(1), -1, IERR )
+            LWORK_ZUNGLQ_M=CDUM(1)
 *           Compute space needed for ZGEBRD
-            CALL ZGEBRD( M, M, A, LDA, S, DUM(1), DUM(1),
-     $                   DUM(1), DUM(1), -1, IERR )
-            LWORK_ZGEBRD=DUM(1)
+            CALL ZGEBRD( M, M, A, LDA, S, DUM(1), CDUM(1),
+     $                   CDUM(1), CDUM(1), -1, IERR )
+            LWORK_ZGEBRD=CDUM(1)
 *            Compute space needed for ZUNGBR P
-            CALL ZUNGBR( 'P', M, M, M, A, N, DUM(1),
-     $                   DUM(1), -1, IERR )
-            LWORK_ZUNGBR_P=DUM(1)
+            CALL ZUNGBR( 'P', M, M, M, A, N, CDUM(1),
+     $                   CDUM(1), -1, IERR )
+            LWORK_ZUNGBR_P=CDUM(1)
 *           Compute space needed for ZUNGBR Q
-            CALL ZUNGBR( 'Q', M, M, M, A, N, DUM(1),
-     $                   DUM(1), -1, IERR )
-            LWORK_ZUNGBR_Q=DUM(1)
+            CALL ZUNGBR( 'Q', M, M, M, A, N, CDUM(1),
+     $                   CDUM(1), -1, IERR )
+            LWORK_ZUNGBR_Q=CDUM(1)
             IF( N.GE.MNTHR ) THEN
                IF( WNTVN ) THEN
 *
@@ -592,21 +593,21 @@
 *
 *              Path 10t(N greater than M, but not much larger)
 *
-               CALL ZGEBRD( M, N, A, LDA, S, DUM(1), DUM(1),
-     $                   DUM(1), DUM(1), -1, IERR )
-               LWORK_ZGEBRD=DUM(1)
+               CALL ZGEBRD( M, N, A, LDA, S, DUM(1), CDUM(1),
+     $                   CDUM(1), CDUM(1), -1, IERR )
+               LWORK_ZGEBRD=CDUM(1)
                MAXWRK = 2*M + LWORK_ZGEBRD
                IF( WNTVS .OR. WNTVO ) THEN
 *                Compute space needed for ZUNGBR P
-                 CALL ZUNGBR( 'P', M, N, M, A, N, DUM(1),
-     $                   DUM(1), -1, IERR )
-                 LWORK_ZUNGBR_P=DUM(1)
+                 CALL ZUNGBR( 'P', M, N, M, A, N, CDUM(1),
+     $                   CDUM(1), -1, IERR )
+                 LWORK_ZUNGBR_P=CDUM(1)
                  MAXWRK = MAX( MAXWRK, 2*M+LWORK_ZUNGBR_P )
                END IF
                IF( WNTVA ) THEN
-                 CALL ZUNGBR( 'P', N,  N, M, A, N, DUM(1),
-     $                   DUM(1), -1, IERR )
-                 LWORK_ZUNGBR_P=DUM(1)
+                 CALL ZUNGBR( 'P', N,  N, M, A, N, CDUM(1),
+     $                   CDUM(1), -1, IERR )
+                 LWORK_ZUNGBR_P=CDUM(1)
                  MAXWRK = MAX( MAXWRK, 2*M+LWORK_ZUNGBR_P )
                END IF
                IF( .NOT.WNTUN ) THEN