diff --git a/SRC/claed0.f b/SRC/claed0.f index f13e000a4..4a35a1bb9 100644 --- a/SRC/claed0.f +++ b/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 * diff --git a/SRC/zlaed0.f b/SRC/zlaed0.f index 3a725dfc4..710c5b75c 100644 --- a/SRC/zlaed0.f +++ b/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 * diff --git a/TESTING/EIG/cchkst.f b/TESTING/EIG/cchkst.f index 790c7f5f1..912e361ba 100644 --- a/TESTING/EIG/cchkst.f +++ b/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, ')' ) * diff --git a/TESTING/EIG/zchkst.f b/TESTING/EIG/zchkst.f index 4335e15f0..dd3ea7a0c 100644 --- a/TESTING/EIG/zchkst.f +++ b/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, ')' ) *