diff --git a/SRC/zgesvj.f b/SRC/zgesvj.f index 311ddf541..85765067f 100644 --- a/SRC/zgesvj.f +++ b/SRC/zgesvj.f @@ -474,7 +474,7 @@ RETURN ELSE IF( LQUERY ) THEN CWORK( 1 ) = LWMIN - RWORK( 1 ) = REAL( LRWMIN ) + RWORK( 1 ) = DBLE( LRWMIN ) RETURN END IF * diff --git a/SRC/zhbevd.f b/SRC/zhbevd.f index 73248d78c..94747ac3a 100644 --- a/SRC/zhbevd.f +++ b/SRC/zhbevd.f @@ -290,7 +290,7 @@ * IF( INFO.EQ.0 ) THEN WORK( 1 ) = LWMIN - RWORK( 1 ) = REAL( LRWMIN ) + RWORK( 1 ) = DBLE( LRWMIN ) IWORK( 1 ) = LIWMIN * IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN @@ -387,7 +387,7 @@ END IF * WORK( 1 ) = LWMIN - RWORK( 1 ) = REAL( LRWMIN ) + RWORK( 1 ) = DBLE( LRWMIN ) IWORK( 1 ) = LIWMIN RETURN * diff --git a/SRC/zhbevd_2stage.f b/SRC/zhbevd_2stage.f index f7036aedd..eda9cf872 100644 --- a/SRC/zhbevd_2stage.f +++ b/SRC/zhbevd_2stage.f @@ -344,7 +344,7 @@ * IF( INFO.EQ.0 ) THEN WORK( 1 ) = LWMIN - RWORK( 1 ) = REAL( LRWMIN ) + RWORK( 1 ) = DBLE( LRWMIN ) IWORK( 1 ) = LIWMIN * IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN @@ -446,7 +446,7 @@ END IF * WORK( 1 ) = LWMIN - RWORK( 1 ) = REAL( LRWMIN ) + RWORK( 1 ) = DBLE( LRWMIN ) IWORK( 1 ) = LIWMIN RETURN * diff --git a/SRC/zhbgvd.f b/SRC/zhbgvd.f index 6b3279faa..82d62622c 100644 --- a/SRC/zhbgvd.f +++ b/SRC/zhbgvd.f @@ -324,7 +324,7 @@ * IF( INFO.EQ.0 ) THEN WORK( 1 ) = LWMIN - RWORK( 1 ) = REAL( LRWMIN ) + RWORK( 1 ) = DBLE( LRWMIN ) IWORK( 1 ) = LIWMIN * IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN @@ -391,7 +391,7 @@ END IF * WORK( 1 ) = LWMIN - RWORK( 1 ) = REAL( LRWMIN ) + RWORK( 1 ) = DBLE( LRWMIN ) IWORK( 1 ) = LIWMIN RETURN * diff --git a/SRC/zheevd.f b/SRC/zheevd.f index f8f69417e..01ad3b25c 100644 --- a/SRC/zheevd.f +++ b/SRC/zheevd.f @@ -284,7 +284,7 @@ LIOPT = LIWMIN END IF WORK( 1 ) = LOPT - RWORK( 1 ) = REAL( LROPT ) + RWORK( 1 ) = DBLE( LROPT ) IWORK( 1 ) = LIOPT * IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN @@ -380,7 +380,7 @@ END IF * WORK( 1 ) = LOPT - RWORK( 1 ) = REAL( LROPT ) + RWORK( 1 ) = DBLE( LROPT ) IWORK( 1 ) = LIOPT * RETURN diff --git a/SRC/zheevd_2stage.f b/SRC/zheevd_2stage.f index 216109914..b73d7155e 100644 --- a/SRC/zheevd_2stage.f +++ b/SRC/zheevd_2stage.f @@ -337,7 +337,7 @@ END IF END IF WORK( 1 ) = LWMIN - RWORK( 1 ) = REAL( LRWMIN ) + RWORK( 1 ) = DBLE( LRWMIN ) IWORK( 1 ) = LIWMIN * IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN @@ -436,7 +436,7 @@ END IF * WORK( 1 ) = LWMIN - RWORK( 1 ) = REAL( LRWMIN ) + RWORK( 1 ) = DBLE( LRWMIN ) IWORK( 1 ) = LIWMIN * RETURN diff --git a/SRC/zheevr.f b/SRC/zheevr.f index 603424126..038738ec8 100644 --- a/SRC/zheevr.f +++ b/SRC/zheevr.f @@ -479,7 +479,7 @@ NB = MAX( NB, ILAENV( 1, 'ZUNMTR', UPLO, N, -1, -1, -1 ) ) LWKOPT = MAX( ( NB+1 )*N, LWMIN ) WORK( 1 ) = LWKOPT - RWORK( 1 ) = REAL( LRWMIN ) + RWORK( 1 ) = DBLE( LRWMIN ) IWORK( 1 ) = LIWMIN * IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN @@ -736,7 +736,7 @@ * Set WORK(1) to optimal workspace size. * WORK( 1 ) = LWKOPT - RWORK( 1 ) = REAL( LRWMIN ) + RWORK( 1 ) = DBLE( LRWMIN ) IWORK( 1 ) = LIWMIN * RETURN diff --git a/SRC/zheevr_2stage.f b/SRC/zheevr_2stage.f index b520f1228..0ba5a2953 100644 --- a/SRC/zheevr_2stage.f +++ b/SRC/zheevr_2stage.f @@ -521,7 +521,7 @@ * IF( INFO.EQ.0 ) THEN WORK( 1 ) = LWMIN - RWORK( 1 ) = REAL( LRWMIN ) + RWORK( 1 ) = DBLE( LRWMIN ) IWORK( 1 ) = LIWMIN * IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN @@ -781,7 +781,7 @@ * Set WORK(1) to optimal workspace size. * WORK( 1 ) = LWMIN - RWORK( 1 ) = REAL( LRWMIN ) + RWORK( 1 ) = DBLE( LRWMIN ) IWORK( 1 ) = LIWMIN * RETURN diff --git a/SRC/zhpevd.f b/SRC/zhpevd.f index 540469a57..1a033de79 100644 --- a/SRC/zhpevd.f +++ b/SRC/zhpevd.f @@ -270,7 +270,7 @@ END IF END IF WORK( 1 ) = LWMIN - RWORK( 1 ) = REAL( LRWMIN ) + RWORK( 1 ) = DBLE( LRWMIN ) IWORK( 1 ) = LIWMIN * IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN @@ -363,7 +363,7 @@ END IF * WORK( 1 ) = LWMIN - RWORK( 1 ) = REAL( LRWMIN ) + RWORK( 1 ) = DBLE( LRWMIN ) IWORK( 1 ) = LIWMIN RETURN * diff --git a/SRC/zstedc.f b/SRC/zstedc.f index 16d90a216..4a5d9fa69 100644 --- a/SRC/zstedc.f +++ b/SRC/zstedc.f @@ -296,7 +296,7 @@ LIWMIN = 3 + 5*N END IF WORK( 1 ) = LWMIN - RWORK( 1 ) = REAL( LRWMIN ) + RWORK( 1 ) = DBLE( LRWMIN ) IWORK( 1 ) = LIWMIN * IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN @@ -472,7 +472,7 @@ * 70 CONTINUE WORK( 1 ) = LWMIN - RWORK( 1 ) = REAL( LRWMIN ) + RWORK( 1 ) = DBLE( LRWMIN ) IWORK( 1 ) = LIWMIN * RETURN diff --git a/TESTING/CMakeLists.txt b/TESTING/CMakeLists.txt index f9c2482fc..58164e516 100644 --- a/TESTING/CMakeLists.txt +++ b/TESTING/CMakeLists.txt @@ -177,6 +177,9 @@ add_lapack_test(zlse.out lse.in xeigtstz) # # ======== COMPLEX16 DMD EIG TESTS =========================== add_lapack_test(zdmd.out zdmd.in xdmdeigtstz) +# +# ======== COMPLEX16 WORKSPACE QUERY PRECISION TEST =========== +add_test(NAME LAPACK-test_wq_zrwork COMMAND $) endif() diff --git a/TESTING/EIG/CMakeLists.txt b/TESTING/EIG/CMakeLists.txt index d99762d43..66ead8798 100644 --- a/TESTING/EIG/CMakeLists.txt +++ b/TESTING/EIG/CMakeLists.txt @@ -128,4 +128,5 @@ endif() if(BUILD_COMPLEX16) add_eig_executable(xeigtstz ${ZEIGTST} ${DZIGTST} ${AEIGTST}) add_eig_executable(xdmdeigtstz ${ZDMDEIGTST}) +add_eig_executable(test_wq_zrwork test_wq_rwork.f) endif() diff --git a/TESTING/EIG/test_wq_rwork.f b/TESTING/EIG/test_wq_rwork.f new file mode 100644 index 000000000..efaa0701e --- /dev/null +++ b/TESTING/EIG/test_wq_rwork.f @@ -0,0 +1,182 @@ +*> \brief Test workspace query precision for z* RWORK +* +* =========== DOCUMENTATION =========== +* +* Purpose +* ======= +* +* TEST_WQ_RWORK validates that workspace query calls (LWORK=-1, +* LRWORK=-1, LIWORK=-1) return exact RWORK sizes for COMPLEX*16 +* routines with O(N^2) LRWMIN formulas. +* +* When LRWMIN > 2^24 (approx N > 2896 for formula 1+5N+2N^2), +* storing the value through a REAL (float32) intermediary loses +* precision. This test catches that regression by checking +* INT(RWORK(1)) == expected at N values above the threshold. +* +* No large matrices are allocated -- workspace queries return +* immediately after storing sizes, so the test runs in microseconds. +* +* =========== END DOCUMENTATION ======== +* + PROGRAM TEST_WQ_RWORK +* + IMPLICIT NONE +* +* .. Parameters .. + INTEGER NNVALS + PARAMETER ( NNVALS = 3 ) +* +* .. Local Scalars .. + INTEGER INFO, N, LRWEXP, NFAIL, NPASS, I, LRWGOT +* +* .. Local Arrays .. +* Minimal dummy arrays for workspace queries (never accessed +* by the routines when LWORK=-1). + INTEGER NVALS( NNVALS ), IWORK( 1 ) + COMPLEX*16 A( 1 ), B( 1 ), AB( 1 ), BB( 1 ) + COMPLEX*16 AP( 1 ), BP( 1 ), Z( 1 ), WORK( 1 ) + DOUBLE PRECISION W( 1 ), RWORK( 1 ), D( 1 ), E( 1 ) +* +* .. External Subroutines .. + EXTERNAL ZHEEVD, ZHEGVD, ZHBEVD, ZHPEVD + EXTERNAL ZHPGVD, ZHBGVD, ZSTEDC +* +* Test N values: 1000 (below 2^24 threshold), 3000 and 5000 (above) + DATA NVALS / 1000, 3000, 5000 / +* +* .. Executable Statements .. +* + NFAIL = 0 + NPASS = 0 +* + WRITE( *, * ) 'Workspace query precision test for z* RWORK' + WRITE( *, * ) '============================================' + WRITE( *, * ) +* + DO 100 I = 1, NNVALS + N = NVALS( I ) +* +* Expected LRWMIN for JOBZ='V': 1 + 5*N + 2*N**2 +* (common to ZHEEVD, ZHEGVD, ZHBEVD, ZHPEVD, ZHPGVD, ZHBGVD) +* + LRWEXP = 1 + 5*N + 2*N*N +* +* ---- ZHEEVD ---- +* + INFO = 0 + CALL ZHEEVD( 'V', 'U', N, A, N, W, + $ WORK, -1, RWORK, -1, IWORK, -1, INFO ) + LRWGOT = INT( RWORK( 1 ) ) + IF( INFO.EQ.0 .AND. LRWGOT.EQ.LRWEXP ) THEN + NPASS = NPASS + 1 + ELSE + NFAIL = NFAIL + 1 + WRITE( *, 9999 ) 'ZHEEVD ', N, LRWEXP, LRWGOT, INFO + END IF +* +* ---- ZHEGVD ---- +* + INFO = 0 + CALL ZHEGVD( 1, 'V', 'U', N, A, N, B, N, W, + $ WORK, -1, RWORK, -1, IWORK, -1, INFO ) + LRWGOT = INT( RWORK( 1 ) ) + IF( INFO.EQ.0 .AND. LRWGOT.EQ.LRWEXP ) THEN + NPASS = NPASS + 1 + ELSE + NFAIL = NFAIL + 1 + WRITE( *, 9999 ) 'ZHEGVD ', N, LRWEXP, LRWGOT, INFO + END IF +* +* ---- ZHBEVD ---- +* KD=0 (diagonal band matrix), LDAB=1, LDZ=N +* + INFO = 0 + CALL ZHBEVD( 'V', 'U', N, 0, AB, 1, W, Z, N, + $ WORK, -1, RWORK, -1, IWORK, -1, INFO ) + LRWGOT = INT( RWORK( 1 ) ) + IF( INFO.EQ.0 .AND. LRWGOT.EQ.LRWEXP ) THEN + NPASS = NPASS + 1 + ELSE + NFAIL = NFAIL + 1 + WRITE( *, 9999 ) 'ZHBEVD ', N, LRWEXP, LRWGOT, INFO + END IF +* +* ---- ZHPEVD ---- +* LDZ=N +* + INFO = 0 + CALL ZHPEVD( 'V', 'U', N, AP, W, Z, N, + $ WORK, -1, RWORK, -1, IWORK, -1, INFO ) + LRWGOT = INT( RWORK( 1 ) ) + IF( INFO.EQ.0 .AND. LRWGOT.EQ.LRWEXP ) THEN + NPASS = NPASS + 1 + ELSE + NFAIL = NFAIL + 1 + WRITE( *, 9999 ) 'ZHPEVD ', N, LRWEXP, LRWGOT, INFO + END IF +* +* ---- ZHPGVD ---- +* ITYPE=1, LDZ=N +* + INFO = 0 + CALL ZHPGVD( 1, 'V', 'U', N, AP, BP, W, Z, N, + $ WORK, -1, RWORK, -1, IWORK, -1, INFO ) + LRWGOT = INT( RWORK( 1 ) ) + IF( INFO.EQ.0 .AND. LRWGOT.EQ.LRWEXP ) THEN + NPASS = NPASS + 1 + ELSE + NFAIL = NFAIL + 1 + WRITE( *, 9999 ) 'ZHPGVD ', N, LRWEXP, LRWGOT, INFO + END IF +* +* ---- ZHBGVD ---- +* KA=0, KB=0, LDAB=1, LDBB=1, LDZ=N +* + INFO = 0 + CALL ZHBGVD( 'V', 'U', N, 0, 0, AB, 1, BB, 1, + $ W, Z, N, + $ WORK, -1, RWORK, -1, IWORK, -1, INFO ) + LRWGOT = INT( RWORK( 1 ) ) + IF( INFO.EQ.0 .AND. LRWGOT.EQ.LRWEXP ) THEN + NPASS = NPASS + 1 + ELSE + NFAIL = NFAIL + 1 + WRITE( *, 9999 ) 'ZHBGVD ', N, LRWEXP, LRWGOT, INFO + END IF +* +* ---- ZSTEDC (COMPZ='I') ---- +* Expected: 1 + 4*N + 2*N**2 +* + LRWEXP = 1 + 4*N + 2*N*N +* + INFO = 0 + CALL ZSTEDC( 'I', N, D, E, Z, N, + $ WORK, -1, RWORK, -1, IWORK, -1, INFO ) + LRWGOT = INT( RWORK( 1 ) ) + IF( INFO.EQ.0 .AND. LRWGOT.EQ.LRWEXP ) THEN + NPASS = NPASS + 1 + ELSE + NFAIL = NFAIL + 1 + WRITE( *, 9999 ) 'ZSTEDC ', N, LRWEXP, LRWGOT, INFO + END IF +* + 100 CONTINUE +* +* Print summary +* + WRITE( *, * ) + IF( NFAIL.EQ.0 ) THEN + WRITE( *, 9998 ) NPASS + ELSE + WRITE( *, 9997 ) NFAIL, NFAIL + NPASS + END IF +* + 9999 FORMAT( ' FAIL: ', A7, ' N=', I6, + $ ' expected=', I12, ' got=', I12, ' INFO=', I4 ) + 9998 FORMAT( ' All ', I3, ' workspace query tests PASSED' ) + 9997 FORMAT( ' ', I3, ' of ', I3, ' tests FAILED' ) +* + IF( NFAIL.NE.0 ) STOP 1 +* + END