Keep the divide and conquer workspace of cLAED0 and zLAED0 out of the caller's Z

SRC/claed0.f

@@ -105,16 +105,15 @@
 *> \param[out] QSTORE
 *> \verbatim
 *>          QSTORE is COMPLEX array, dimension (LDQS, N)
-*>         Used to store parts of
-*>         the eigenvector matrix when the updating matrix multiplies
-*>         take place.
+*>         Workspace: holds the products of the leaf eigenvector
+*>         matrices and the QSIZ*N workspace of CLAED7.
 *> \endverbatim
 *>
 *> \param[in] LDQS
 *> \verbatim
 *>          LDQS is INTEGER
 *>         The leading dimension of the array QSTORE.
-*>         LDQS >= max(1,N).
+*>         LDQS >= max(1,QSIZ).
 *> \endverbatim
 *>
 *> \param[out] INFO
@@ -171,7 +170,7 @@ SUBROUTINE CLAED0( QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK,
       REAL               TEMP
 *     ..
 *     .. External Subroutines ..
-      EXTERNAL           CCOPY, CLACRM, CLAED7, SCOPY, SSTEQR,
+      EXTERNAL           CCOPY, CLACPY, CLACRM, CLAED7, SCOPY, SSTEQR,
      $                   XERBLA
 *     ..
 *     .. External Functions ..
@@ -288,6 +287,8 @@ SUBROUTINE CLAED0( QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK,
          CALL CLACRM( QSIZ, MATSIZ, Q( 1, SUBMAT ), LDQ, RWORK( LL ),
      $                MATSIZ, QSTORE( 1, SUBMAT ), LDQS,
      $                RWORK( IWREM ) )
+         CALL CLACPY( 'A', QSIZ, MATSIZ, QSTORE( 1, SUBMAT ), LDQS,
+     $                Q( 1, SUBMAT ), LDQ )
          IWORK( IQPTR+CURR+1 ) = IWORK( IQPTR+CURR ) + MATSIZ**2
          CURR = CURR + 1
          IF( INFO.GT.0 ) THEN
@@ -328,15 +329,17 @@ SUBROUTINE CLAED0( QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK,
 *     when the eigenvectors of a full or band Hermitian matrix (which
 *     was reduced to tridiagonal form) are desired.
 *
-*     I am free to use Q as a valuable working space until Loop 150.
+*     The eigenvectors stay in Q, which has leading dimension LDQ, and
+*     the contiguous QSTORE is the QSIZ*MATSIZ workspace of CLAED7:
+*     Q(N+1:LDQ,:) belongs to the caller and must not be touched.
 *
             CALL CLAED7( MATSIZ, MSD2, QSIZ, TLVLS, CURLVL, CURPRB,
-     $                   D( SUBMAT ), QSTORE( 1, SUBMAT ), LDQS,
+     $                   D( SUBMAT ), Q( 1, SUBMAT ), LDQ,
      $                   E( SUBMAT+MSD2-1 ), IWORK( INDXQ+SUBMAT ),
      $                   RWORK( IQ ), IWORK( IQPTR ), IWORK( IPRMPT ),
      $                   IWORK( IPERM ), IWORK( IGIVPT ),
      $                   IWORK( IGIVCL ), RWORK( IGIVNM ),
-     $                   Q( 1, SUBMAT ), RWORK( IWREM ),
+     $                   QSTORE, RWORK( IWREM ),
      $                   IWORK( SUBPBS+1 ), INFO )
             IF( INFO.GT.0 ) THEN
                INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1
@@ -357,9 +360,10 @@ SUBROUTINE CLAED0( QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK,
       DO 100 I = 1, N
          J = IWORK( INDXQ+I )
          RWORK( I ) = D( J )
-         CALL CCOPY( QSIZ, QSTORE( 1, J ), 1, Q( 1, I ), 1 )
+         CALL CCOPY( QSIZ, Q( 1, J ), 1, QSTORE( 1, I ), 1 )
   100 CONTINUE
       CALL SCOPY( N, RWORK, 1, D, 1 )
+      CALL CLACPY( 'A', QSIZ, N, QSTORE, LDQS, Q, LDQ )
 *
       RETURN
 *

SRC/zlaed0.f

@@ -105,16 +105,15 @@
 *> \param[out] QSTORE
 *> \verbatim
 *>          QSTORE is COMPLEX*16 array, dimension (LDQS, N)
-*>         Used to store parts of
-*>         the eigenvector matrix when the updating matrix multiplies
-*>         take place.
+*>         Workspace: holds the products of the leaf eigenvector
+*>         matrices and the QSIZ*N workspace of ZLAED7.
 *> \endverbatim
 *>
 *> \param[in] LDQS
 *> \verbatim
 *>          LDQS is INTEGER
 *>         The leading dimension of the array QSTORE.
-*>         LDQS >= max(1,N).
+*>         LDQS >= max(1,QSIZ).
 *> \endverbatim
 *>
 *> \param[out] INFO
@@ -171,7 +170,7 @@ SUBROUTINE ZLAED0( QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK,
       DOUBLE PRECISION   TEMP
 *     ..
 *     .. External Subroutines ..
-      EXTERNAL           DCOPY, DSTEQR, XERBLA, ZCOPY, ZLACRM,
+      EXTERNAL           DCOPY, DSTEQR, XERBLA, ZCOPY, ZLACPY, ZLACRM,
      $                   ZLAED7
 *     ..
 *     .. External Functions ..
@@ -288,6 +287,8 @@ SUBROUTINE ZLAED0( QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK,
          CALL ZLACRM( QSIZ, MATSIZ, Q( 1, SUBMAT ), LDQ, RWORK( LL ),
      $                MATSIZ, QSTORE( 1, SUBMAT ), LDQS,
      $                RWORK( IWREM ) )
+         CALL ZLACPY( 'A', QSIZ, MATSIZ, QSTORE( 1, SUBMAT ), LDQS,
+     $                Q( 1, SUBMAT ), LDQ )
          IWORK( IQPTR+CURR+1 ) = IWORK( IQPTR+CURR ) + MATSIZ**2
          CURR = CURR + 1
          IF( INFO.GT.0 ) THEN
@@ -328,15 +329,17 @@ SUBROUTINE ZLAED0( QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK,
 *     when the eigenvectors of a full or band Hermitian matrix (which
 *     was reduced to tridiagonal form) are desired.
 *
-*     I am free to use Q as a valuable working space until Loop 150.
+*     The eigenvectors stay in Q, which has leading dimension LDQ, and
+*     the contiguous QSTORE is the QSIZ*MATSIZ workspace of ZLAED7:
+*     Q(N+1:LDQ,:) belongs to the caller and must not be touched.
 *
             CALL ZLAED7( MATSIZ, MSD2, QSIZ, TLVLS, CURLVL, CURPRB,
-     $                   D( SUBMAT ), QSTORE( 1, SUBMAT ), LDQS,
+     $                   D( SUBMAT ), Q( 1, SUBMAT ), LDQ,
      $                   E( SUBMAT+MSD2-1 ), IWORK( INDXQ+SUBMAT ),
      $                   RWORK( IQ ), IWORK( IQPTR ), IWORK( IPRMPT ),
      $                   IWORK( IPERM ), IWORK( IGIVPT ),
      $                   IWORK( IGIVCL ), RWORK( IGIVNM ),
-     $                   Q( 1, SUBMAT ), RWORK( IWREM ),
+     $                   QSTORE, RWORK( IWREM ),
      $                   IWORK( SUBPBS+1 ), INFO )
             IF( INFO.GT.0 ) THEN
                INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1
@@ -357,9 +360,10 @@ SUBROUTINE ZLAED0( QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK,
       DO 100 I = 1, N
          J = IWORK( INDXQ+I )
          RWORK( I ) = D( J )
-         CALL ZCOPY( QSIZ, QSTORE( 1, J ), 1, Q( 1, I ), 1 )
+         CALL ZCOPY( QSIZ, Q( 1, J ), 1, QSTORE( 1, I ), 1 )
   100 CONTINUE
       CALL DCOPY( N, RWORK, 1, D, 1 )
+      CALL ZLACPY( 'A', QSIZ, N, QSTORE, LDQS, Q, LDQ )
 *
       RETURN
 *

TESTING/EIG/cchkst.f

@@ -651,11 +651,14 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
      $                   RTUNFL, TEMP1, TEMP2, TEMP3, TEMP4, ULP,
      $                   ULPINV, UNFL, VL, VU
 *     ..
+      INTEGER            LDZREG, NSMLSZ
 *     .. Local Arrays ..
       INTEGER            IDUMMA( 1 ), IOLDSD( 4 ), ISEED2( 4 ),
      $                   KMAGN( MAXTYP ), KMODE( MAXTYP ),
      $                   KTYPE( MAXTYP )
       REAL               DUMMA( 1 )
+      REAL, ALLOCATABLE :: DREG( : ), EREG( : ), D1REG( : )
+      COMPLEX, ALLOCATABLE :: ZREG( :, : )
 *     ..
 *     .. External Functions ..
       INTEGER            ILAENV
@@ -1953,11 +1956,76 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
   300    CONTINUE
   310 CONTINUE
 *
+*
+*     CSTEDC( 'V' ) must leave the rows of Z below N alone.  CLAED0
+*     used the caller's Z as a packed QSIZ by N workspace, which
+*     spills into rows N+1 to LDZ whenever LDZ > N.  The sizes in the
+*     input file do not reach the divide and conquer recursion, so it
+*     is exercised with local arrays and N > SMLSIZ.  LDU does not
+*     give the capacity of the caller's eigenvalue arrays or Z columns.
+*
+      NSMLSZ = ILAENV( 9, 'CSTEDC', ' ', 0, 0, 0, 0 )
+      N = MAX( 2, NSMLSZ+1 )
+      LDZREG = N + 1
+      ALLOCATE( DREG( N ), EREG( N ), D1REG( N ),
+     $          ZREG( LDZREG, N ) )
+      CALL CSTEDC( 'V', N, DREG, RWORK, ZREG, LDZREG, WORK, -1,
+     $             RWORK, -1, IWORK, -1, IINFO )
+      IF( IINFO.EQ.0 .AND. INT( REAL( WORK( 1 ) ) ).LE.LWORK .AND.
+     $    INT( RWORK( 1 ) ).LE.LRWORK-N .AND. IWORK( 1 ).LE.LIWORK )
+     $    THEN
+         DO 390 J = 1, N
+            DREG( J ) = REAL( J )
+            EREG( J ) = ONE / REAL( J+1 )
+  390    CONTINUE
+         EREG( N ) = ZERO
+         CALL SCOPY( N, DREG, 1, D1REG, 1 )
+         CALL SCOPY( N-1, EREG, 1, RWORK, 1 )
+         CALL CLASET( 'Full', N, N, CZERO, CONE, ZREG, LDZREG )
+         CALL CLASET( 'Full', LDZREG-N, N, -CONE, -CONE,
+     $                ZREG( N+1, 1 ), LDZREG )
+         CALL CSTEDC( 'V', N, D1REG, RWORK, ZREG, LDZREG, WORK, LWORK,
+     $                RWORK( N+1 ), LRWORK-N, IWORK, LIWORK, IINFO )
+         NTESTT = NTESTT + 1
+         IF( IINFO.NE.0 ) THEN
+            WRITE( NOUNIT, FMT = 9982 )N, IINFO
+            NERRS = NERRS + 1
+         ELSE
+            ITEMP = 0
+            DO 410 J = 1, N
+               DO 400 I = N + 1, LDZREG
+                  IF( ZREG( I, J ).NE.-CONE )
+     $               ITEMP = ITEMP + 1
+  400          CONTINUE
+  410       CONTINUE
+            IF( ITEMP.GT.0 ) THEN
+               WRITE( NOUNIT, FMT = 9981 )N, ITEMP
+               NERRS = NERRS + 1
+            END IF
+            CALL CSTT21( N, 0, DREG, EREG, D1REG, DUMMA, ZREG,
+     $                   LDZREG, WORK, RWORK( N+1 ), RESULT( 1 ) )
+            DO 420 J = 1, 2
+               IF( RESULT( J ).GE.THRESH ) THEN
+                  WRITE( NOUNIT, FMT = 9980 )N, J, RESULT( J )
+                  NERRS = NERRS + 1
+               END IF
+  420       CONTINUE
+            NTESTT = NTESTT + 3
+         END IF
+      END IF
+      DEALLOCATE( DREG, EREG, D1REG, ZREG )
+*
 *     Summary
 *
       CALL SLASUM( 'CST', NOUNIT, NERRS, NTESTT )
       RETURN
 *
+ 9982 FORMAT( ' CCHKST: CSTEDC( V ) with N=', I5, ' returned INFO=',
+     $      I6 )
+ 9981 FORMAT( ' CCHKST: CSTEDC( V ) with N=', I5, ' overwrote ', I6,
+     $      ' entries of Z below row N' )
+ 9980 FORMAT( ' CCHKST: CSTEDC( V ) with N=', I5, ', test ', I2,
+     $      ' ratio=', G10.3, ' >= threshold' )
  9999 FORMAT( ' CCHKST: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
      $      I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
 *

TESTING/EIG/zchkst.f

@@ -651,11 +651,14 @@ SUBROUTINE ZCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
      $                   RTUNFL, TEMP1, TEMP2, TEMP3, TEMP4, ULP,
      $                   ULPINV, UNFL, VL, VU
 *     ..
+      INTEGER            LDZREG, NSMLSZ
 *     .. Local Arrays ..
       INTEGER            IDUMMA( 1 ), IOLDSD( 4 ), ISEED2( 4 ),
      $                   KMAGN( MAXTYP ), KMODE( MAXTYP ),
      $                   KTYPE( MAXTYP )
       DOUBLE PRECISION   DUMMA( 1 )
+      DOUBLE PRECISION, ALLOCATABLE :: DREG( : ), EREG( : ), D1REG( : )
+      COMPLEX*16, ALLOCATABLE :: ZREG( :, : )
 *     ..
 *     .. External Functions ..
       INTEGER            ILAENV
@@ -1952,11 +1955,76 @@ SUBROUTINE ZCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
   300    CONTINUE
   310 CONTINUE
 *
+*
+*     ZSTEDC( 'V' ) must leave the rows of Z below N alone.  ZLAED0
+*     used the caller's Z as a packed QSIZ by N workspace, which
+*     spills into rows N+1 to LDZ whenever LDZ > N.  The sizes in the
+*     input file do not reach the divide and conquer recursion, so it
+*     is exercised with local arrays and N > SMLSIZ.  LDU does not
+*     give the capacity of the caller's eigenvalue arrays or Z columns.
+*
+      NSMLSZ = ILAENV( 9, 'ZSTEDC', ' ', 0, 0, 0, 0 )
+      N = MAX( 2, NSMLSZ+1 )
+      LDZREG = N + 1
+      ALLOCATE( DREG( N ), EREG( N ), D1REG( N ),
+     $          ZREG( LDZREG, N ) )
+      CALL ZSTEDC( 'V', N, DREG, RWORK, ZREG, LDZREG, WORK, -1,
+     $             RWORK, -1, IWORK, -1, IINFO )
+      IF( IINFO.EQ.0 .AND. INT( DBLE( WORK( 1 ) ) ).LE.LWORK .AND.
+     $    INT( RWORK( 1 ) ).LE.LRWORK-N .AND. IWORK( 1 ).LE.LIWORK )
+     $    THEN
+         DO 390 J = 1, N
+            DREG( J ) = DBLE( J )
+            EREG( J ) = ONE / DBLE( J+1 )
+  390    CONTINUE
+         EREG( N ) = ZERO
+         CALL DCOPY( N, DREG, 1, D1REG, 1 )
+         CALL DCOPY( N-1, EREG, 1, RWORK, 1 )
+         CALL ZLASET( 'Full', N, N, CZERO, CONE, ZREG, LDZREG )
+         CALL ZLASET( 'Full', LDZREG-N, N, -CONE, -CONE,
+     $                ZREG( N+1, 1 ), LDZREG )
+         CALL ZSTEDC( 'V', N, D1REG, RWORK, ZREG, LDZREG, WORK, LWORK,
+     $                RWORK( N+1 ), LRWORK-N, IWORK, LIWORK, IINFO )
+         NTESTT = NTESTT + 1
+         IF( IINFO.NE.0 ) THEN
+            WRITE( NOUNIT, FMT = 9982 )N, IINFO
+            NERRS = NERRS + 1
+         ELSE
+            ITEMP = 0
+            DO 410 J = 1, N
+               DO 400 I = N + 1, LDZREG
+                  IF( ZREG( I, J ).NE.-CONE )
+     $               ITEMP = ITEMP + 1
+  400          CONTINUE
+  410       CONTINUE
+            IF( ITEMP.GT.0 ) THEN
+               WRITE( NOUNIT, FMT = 9981 )N, ITEMP
+               NERRS = NERRS + 1
+            END IF
+            CALL ZSTT21( N, 0, DREG, EREG, D1REG, DUMMA, ZREG,
+     $                   LDZREG, WORK, RWORK( N+1 ), RESULT( 1 ) )
+            DO 420 J = 1, 2
+               IF( RESULT( J ).GE.THRESH ) THEN
+                  WRITE( NOUNIT, FMT = 9980 )N, J, RESULT( J )
+                  NERRS = NERRS + 1
+               END IF
+  420       CONTINUE
+            NTESTT = NTESTT + 3
+         END IF
+      END IF
+      DEALLOCATE( DREG, EREG, D1REG, ZREG )
+*
 *     Summary
 *
       CALL DLASUM( 'ZST', NOUNIT, NERRS, NTESTT )
       RETURN
 *
+ 9982 FORMAT( ' ZCHKST: ZSTEDC( V ) with N=', I5, ' returned INFO=',
+     $      I6 )
+ 9981 FORMAT( ' ZCHKST: ZSTEDC( V ) with N=', I5, ' overwrote ', I6,
+     $      ' entries of Z below row N' )
+ 9980 FORMAT( ' ZCHKST: ZSTEDC( V ) with N=', I5, ', test ', I2,
+     $      ' ratio=', G10.3, ' >= threshold' )
  9999 FORMAT( ' ZCHKST: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
      $      I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
 *