From 5b198e089acf112af8f652db8721981f1e31ce4f Mon Sep 17 00:00:00 2001 From: Rasmus Munk Larsen Date: Thu, 10 Sep 2026 01:11:18 -0700 Subject: [PATCH 1/2] Keep the divide and conquer workspace of cLAED0 and zLAED0 out of the caller's Z cLAED0 and zLAED0 solve the leaf problems into QSTORE, the N by N workspace with leading dimension LDQS, and then merge them with xLAED7 while passing the caller's Q, which has leading dimension LDQ, as the QSIZ*N packed workspace of xLAED7: CALL ZLAED7( ..., D( SUBMAT ), QSTORE( 1, SUBMAT ), LDQS, ..., $ Q( 1, SUBMAT ), RWORK( IWREM ), ... ) xLAED7 and xLAED8 address that workspace with leading dimension QSIZ, so whenever LDQ > QSIZ the merge writes into rows QSIZ+1 to LDQ of Q, which belong to the caller. For cSTEDC and zSTEDC with COMPZ = 'V' and N > SMLSIZ that is every entry of Z(N+1:LDZ, :) for a user whose Z is a block of a larger array: with LDZ = N + 8 the rows N+1 to N+8 of the first columns come back overwritten for any matrix. The real routines keep their workspace in WORK and are not affected. Keep the eigenvectors in Q, with its own leading dimension, copying each leaf product back from QSTORE, and hand QSTORE, which is contiguous and holds LDQS*N >= QSIZ*MATSIZ entries, to xLAED7 as its workspace. The final permutation of the columns goes through QSTORE and back into Q. This costs one copy of QSIZ by MATSIZ per leaf and one of QSIZ by N at the end, O(N^2) against the O(N^3) merges, and changes no arithmetic: the eigenvalues and eigenvectors are bit-identical to the parent commit. The documentation of QSTORE and of its leading dimension, LDQS >= max(1,QSIZ), which the leaf products already required, is updated with it. cCHKST and zCHKST get the case as a regression test for the divide and conquer path, which the sizes in sep.in (N <= 20, below SMLSIZ = 25) never reach: after the size and type loops they call xSTEDC with COMPZ = 'V' on an N = LDU - 1 tridiagonal matrix with the row below N of Z set to a marker, and report an overwritten marker, a nonzero INFO, or a residual or orthogonality ratio above the threshold as a failure. On the parent commit the marker row is overwritten in every column in both precisions; here it is intact and the ratios pass. The full LAPACK test suite passes: 0 numerical errors, 0 other errors, 40 tests more than the parent from the new calls. The regression test fails on the parent and passes with the fix, and the reproducer prints the same before and after output, with gfortran 13 (x86-64 Release and Debug with -fcheck=all, and under QEMU on aarch64, ppc64le, s390x and riscv64), flang-19 and Intel ifx 2025.3. Co-Authored-By: Claude Fable 5.1 --- SRC/claed0.f | 22 ++++++++------- SRC/zlaed0.f | 22 ++++++++------- TESTING/EIG/cchkst.f | 64 ++++++++++++++++++++++++++++++++++++++++++++ TESTING/EIG/zchkst.f | 64 ++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 154 insertions(+), 18 deletions(-) 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..0960c7457 100644 --- a/TESTING/EIG/cchkst.f +++ b/TESTING/EIG/cchkst.f @@ -651,6 +651,7 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, $ RTUNFL, TEMP1, TEMP2, TEMP3, TEMP4, ULP, $ ULPINV, UNFL, VL, VU * .. + INTEGER NSMLSZ * .. Local Arrays .. INTEGER IDUMMA( 1 ), IOLDSD( 4 ), ISEED2( 4 ), $ KMAGN( MAXTYP ), KMODE( MAXTYP ), @@ -1953,11 +1954,74 @@ 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 here with N = LDU - 1, and the eigensystem is +* checked as in the size loop. +* + NSMLSZ = ILAENV( 9, 'CSTEDC', ' ', 0, 0, 0, 0 ) + N = LDU - 1 + IF( N.GT.NSMLSZ ) THEN + CALL CSTEDC( 'V', N, SD, RWORK, Z, LDU, 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 + SD( J ) = REAL( J ) + SE( J ) = ONE / REAL( J+1 ) + 390 CONTINUE + SE( N ) = ZERO + CALL SCOPY( N, SD, 1, D1, 1 ) + CALL SCOPY( N-1, SE, 1, RWORK, 1 ) + CALL CLASET( 'Full', N, N, CZERO, CONE, Z, LDU ) + CALL CLASET( 'Full', LDU-N, N, -CONE, -CONE, Z( N+1, 1 ), + $ LDU ) + CALL CSTEDC( 'V', N, D1, RWORK, Z, LDU, 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, LDU + IF( Z( 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, SD, SE, D1, DUMMA, Z, LDU, 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 + END IF +* * 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..d9f231959 100644 --- a/TESTING/EIG/zchkst.f +++ b/TESTING/EIG/zchkst.f @@ -651,6 +651,7 @@ SUBROUTINE ZCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, $ RTUNFL, TEMP1, TEMP2, TEMP3, TEMP4, ULP, $ ULPINV, UNFL, VL, VU * .. + INTEGER NSMLSZ * .. Local Arrays .. INTEGER IDUMMA( 1 ), IOLDSD( 4 ), ISEED2( 4 ), $ KMAGN( MAXTYP ), KMODE( MAXTYP ), @@ -1952,11 +1953,74 @@ 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 here with N = LDU - 1, and the eigensystem is +* checked as in the size loop. +* + NSMLSZ = ILAENV( 9, 'ZSTEDC', ' ', 0, 0, 0, 0 ) + N = LDU - 1 + IF( N.GT.NSMLSZ ) THEN + CALL ZSTEDC( 'V', N, SD, RWORK, Z, LDU, 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 + SD( J ) = DBLE( J ) + SE( J ) = ONE / DBLE( J+1 ) + 390 CONTINUE + SE( N ) = ZERO + CALL DCOPY( N, SD, 1, D1, 1 ) + CALL DCOPY( N-1, SE, 1, RWORK, 1 ) + CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDU ) + CALL ZLASET( 'Full', LDU-N, N, -CONE, -CONE, Z( N+1, 1 ), + $ LDU ) + CALL ZSTEDC( 'V', N, D1, RWORK, Z, LDU, 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, LDU + IF( Z( 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, SD, SE, D1, DUMMA, Z, LDU, 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 + END IF +* * 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, ')' ) * From 973813c59e6565ce2817d01ca30665b04c48c73c Mon Sep 17 00:00:00 2001 From: Rasmus Munk Larsen Date: Thu, 10 Sep 2026 17:51:27 -0700 Subject: [PATCH 2/2] TESTING: Bound the LAED0 padding regression storage Allocate local eigenvalue arrays and a matrix with one padding row, using a problem size above SMLSIZ. LDU is only a row stride and does not establish the caller's array capacities. Preserve padding, residual, and orthogonality checks. Validation: 4 sep/se2 driver runs and 16 AddressSanitizer capacity cases passed. Both precisions retain the expected padding failures against the parent kernels. --- TESTING/EIG/cchkst.f | 94 +++++++++++++++++++++++--------------------- TESTING/EIG/zchkst.f | 94 +++++++++++++++++++++++--------------------- 2 files changed, 98 insertions(+), 90 deletions(-) diff --git a/TESTING/EIG/cchkst.f b/TESTING/EIG/cchkst.f index 0960c7457..912e361ba 100644 --- a/TESTING/EIG/cchkst.f +++ b/TESTING/EIG/cchkst.f @@ -651,12 +651,14 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, $ RTUNFL, TEMP1, TEMP2, TEMP3, TEMP4, ULP, $ ULPINV, UNFL, VL, VU * .. - INTEGER NSMLSZ + 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 @@ -1959,57 +1961,59 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, * 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 here with N = LDU - 1, and the eigensystem is -* checked as in the size loop. +* 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 = LDU - 1 - IF( N.GT.NSMLSZ ) THEN - CALL CSTEDC( 'V', N, SD, RWORK, Z, LDU, 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 - SD( J ) = REAL( J ) - SE( J ) = ONE / REAL( J+1 ) - 390 CONTINUE - SE( N ) = ZERO - CALL SCOPY( N, SD, 1, D1, 1 ) - CALL SCOPY( N-1, SE, 1, RWORK, 1 ) - CALL CLASET( 'Full', N, N, CZERO, CONE, Z, LDU ) - CALL CLASET( 'Full', LDU-N, N, -CONE, -CONE, Z( N+1, 1 ), - $ LDU ) - CALL CSTEDC( 'V', N, D1, RWORK, Z, LDU, WORK, LWORK, - $ RWORK( N+1 ), LRWORK-N, IWORK, LIWORK, IINFO ) - NTESTT = NTESTT + 1 - IF( IINFO.NE.0 ) THEN - WRITE( NOUNIT, FMT = 9982 )N, IINFO + 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 - ELSE - ITEMP = 0 - DO 410 J = 1, N - DO 400 I = N + 1, LDU - IF( Z( I, J ).NE.-CONE ) - $ ITEMP = ITEMP + 1 - 400 CONTINUE - 410 CONTINUE - IF( ITEMP.GT.0 ) THEN - WRITE( NOUNIT, FMT = 9981 )N, ITEMP + 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 - CALL CSTT21( N, 0, SD, SE, D1, DUMMA, Z, LDU, 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 + 420 CONTINUE + NTESTT = NTESTT + 3 END IF END IF + DEALLOCATE( DREG, EREG, D1REG, ZREG ) * * Summary * diff --git a/TESTING/EIG/zchkst.f b/TESTING/EIG/zchkst.f index d9f231959..dd3ea7a0c 100644 --- a/TESTING/EIG/zchkst.f +++ b/TESTING/EIG/zchkst.f @@ -651,12 +651,14 @@ SUBROUTINE ZCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, $ RTUNFL, TEMP1, TEMP2, TEMP3, TEMP4, ULP, $ ULPINV, UNFL, VL, VU * .. - INTEGER NSMLSZ + 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 @@ -1958,57 +1960,59 @@ SUBROUTINE ZCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, * 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 here with N = LDU - 1, and the eigensystem is -* checked as in the size loop. +* 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 = LDU - 1 - IF( N.GT.NSMLSZ ) THEN - CALL ZSTEDC( 'V', N, SD, RWORK, Z, LDU, 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 - SD( J ) = DBLE( J ) - SE( J ) = ONE / DBLE( J+1 ) - 390 CONTINUE - SE( N ) = ZERO - CALL DCOPY( N, SD, 1, D1, 1 ) - CALL DCOPY( N-1, SE, 1, RWORK, 1 ) - CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDU ) - CALL ZLASET( 'Full', LDU-N, N, -CONE, -CONE, Z( N+1, 1 ), - $ LDU ) - CALL ZSTEDC( 'V', N, D1, RWORK, Z, LDU, WORK, LWORK, - $ RWORK( N+1 ), LRWORK-N, IWORK, LIWORK, IINFO ) - NTESTT = NTESTT + 1 - IF( IINFO.NE.0 ) THEN - WRITE( NOUNIT, FMT = 9982 )N, IINFO + 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 - ELSE - ITEMP = 0 - DO 410 J = 1, N - DO 400 I = N + 1, LDU - IF( Z( I, J ).NE.-CONE ) - $ ITEMP = ITEMP + 1 - 400 CONTINUE - 410 CONTINUE - IF( ITEMP.GT.0 ) THEN - WRITE( NOUNIT, FMT = 9981 )N, ITEMP + 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 - CALL ZSTT21( N, 0, SD, SE, D1, DUMMA, Z, LDU, 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 + 420 CONTINUE + NTESTT = NTESTT + 3 END IF END IF + DEALLOCATE( DREG, EREG, D1REG, ZREG ) * * Summary *