Skip to content

Fix truncation of large workspace values in ZHE routines (Reference-LAPACK PR 1195) #5729

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
martin-frbg merged 1 commit into from
Mar 31, 2026
Merged

Fix truncation of large workspace values in ZHE routines (Reference-LAPACK PR 1195) #5729

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
83 changes: 53 additions & 30 deletions lapack-netlib/SRC/zgesvj.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,15 +5,13 @@
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download ZGESVJ + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgesvj.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgesvj.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgesvj.f">
*> [TXT]</a>
*> \endhtmlonly
*
* Definition:
* ===========
Expand Down Expand Up @@ -101,7 +99,7 @@
*> \param[in] M
*> \verbatim
*> M is INTEGER
*> The number of rows of the input matrix A. 1/DLAMCH('E') > M >= 0.
*> The number of rows of the input matrix A. 1/DLAMCH('E') >= M >= 0.
*> \endverbatim
*>
*> \param[in] N
Expand Down Expand Up @@ -217,7 +215,7 @@
*> LWORK >= 1, if MIN(M,N) = 0, and LWORK >= M+N, otherwise.
*>
*> If on entry LWORK = -1, then a workspace query is assumed and
*> no computation is done; CWORK(1) is set to the minial (and optimal)
*> no computation is done; CWORK(1) is set to the minimal (and optimal)
*> length of CWORK.
*> \endverbatim
*>
Expand Down Expand Up @@ -258,7 +256,7 @@
*> LRWORK >= 1, if MIN(M,N) = 0, and LRWORK >= MAX(6,N), otherwise.
*>
*> If on entry LRWORK = -1, then a workspace query is assumed and
*> no computation is done; RWORK(1) is set to the minial (and optimal)
*> no computation is done; RWORK(1) is set to the minimal (and optimal)
*> length of RWORK.
*> \endverbatim
*>
Expand Down Expand Up @@ -414,7 +412,8 @@ SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
* from BLAS
EXTERNAL ZCOPY, ZROT, ZDSCAL, ZSWAP, ZAXPY
* from LAPACK
EXTERNAL DLASCL, ZLASCL, ZLASET, ZLASSQ, XERBLA
EXTERNAL DLASCL, ZLASCL, ZLASET, ZLASSQ,
$ XERBLA
EXTERNAL ZGSVJ0, ZGSVJ1
* ..
* .. Executable Statements ..
Expand All @@ -440,9 +439,13 @@ SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
LQUERY = ( LWORK.EQ.-1 ) .OR. ( LRWORK.EQ.-1 )
IF( .NOT.( UPPER .OR. LOWER .OR. LSAME( JOBA, 'G' ) ) ) THEN
INFO = -1
ELSE IF( .NOT.( LSVEC .OR. UCTOL .OR. LSAME( JOBU, 'N' ) ) ) THEN
ELSE IF( .NOT.( LSVEC .OR.
$ UCTOL .OR.
$ LSAME( JOBU, 'N' ) ) ) THEN
INFO = -2
ELSE IF( .NOT.( RSVEC .OR. APPLV .OR. LSAME( JOBV, 'N' ) ) ) THEN
ELSE IF( .NOT.( RSVEC .OR.
$ APPLV .OR.
$ LSAME( JOBV, 'N' ) ) ) THEN
INFO = -3
ELSE IF( M.LT.0 ) THEN
INFO = -4
Expand All @@ -455,7 +458,7 @@ SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
ELSE IF( ( RSVEC .AND. ( LDV.LT.N ) ) .OR.
$ ( APPLV .AND. ( LDV.LT.MV ) ) ) THEN
INFO = -11
ELSE IF( UCTOL .AND. ( RWORK( 1 ).LE.ONE ) ) THEN
ELSE IF( UCTOL .AND. ( RWORK( 1 ).LT.ONE ) ) THEN
INFO = -12
ELSE IF( LWORK.LT.LWMIN .AND. ( .NOT.LQUERY ) ) THEN
INFO = -13
Expand All @@ -471,7 +474,7 @@ SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
RETURN
ELSE IF( LQUERY ) THEN
CWORK( 1 ) = LWMIN
RWORK( 1 ) = LRWMIN
RWORK( 1 ) = DBLE( LRWMIN )
RETURN
END IF
*
Expand Down Expand Up @@ -785,7 +788,8 @@ SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
$ V( N4*q+1, N4+1 ), LDV, EPSLN, SFMIN, TOL, 1,
$ CWORK( N+1 ), LWORK-N, IERR )
*
CALL ZGSVJ0( JOBV, M, N4, A, LDA, CWORK, SVA, MVL, V, LDV,
CALL ZGSVJ0( JOBV, M, N4, A, LDA, CWORK, SVA, MVL, V,
$ LDV,
$ EPSLN, SFMIN, TOL, 1, CWORK( N+1 ), LWORK-N,
$ IERR )
*
Expand All @@ -797,16 +801,19 @@ SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
ELSE IF( UPPER ) THEN
*
*
CALL ZGSVJ0( JOBV, N4, N4, A, LDA, CWORK, SVA, MVL, V, LDV,
CALL ZGSVJ0( JOBV, N4, N4, A, LDA, CWORK, SVA, MVL, V,
$ LDV,
$ EPSLN, SFMIN, TOL, 2, CWORK( N+1 ), LWORK-N,
$ IERR )
*
CALL ZGSVJ0( JOBV, N2, N4, A( 1, N4+1 ), LDA, CWORK( N4+1 ),
CALL ZGSVJ0( JOBV, N2, N4, A( 1, N4+1 ), LDA,
$ CWORK( N4+1 ),
$ SVA( N4+1 ), MVL, V( N4*q+1, N4+1 ), LDV,
$ EPSLN, SFMIN, TOL, 1, CWORK( N+1 ), LWORK-N,
$ IERR )
*
CALL ZGSVJ1( JOBV, N2, N2, N4, A, LDA, CWORK, SVA, MVL, V,
CALL ZGSVJ1( JOBV, N2, N2, N4, A, LDA, CWORK, SVA, MVL,
$ V,
$ LDV, EPSLN, SFMIN, TOL, 1, CWORK( N+1 ),
$ LWORK-N, IERR )
*
Expand Down Expand Up @@ -960,7 +967,8 @@ SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
T = HALF / THETA
CS = ONE

CALL ZROT( M, A(1,p), 1, A(1,q), 1,
CALL ZROT( M, A(1,p), 1, A(1,q),
$ 1,
$ CS, CONJG(OMPQ)*T )
IF ( RSVEC ) THEN
CALL ZROT( MVL, V(1,p), 1,
Expand Down Expand Up @@ -989,7 +997,8 @@ SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
AAPP = AAPP*SQRT( MAX( ZERO,
$ ONE-T*AQOAP*AAPQ1 ) )
*
CALL ZROT( M, A(1,p), 1, A(1,q), 1,
CALL ZROT( M, A(1,p), 1, A(1,q),
$ 1,
$ CS, CONJG(OMPQ)*SN )
IF ( RSVEC ) THEN
CALL ZROT( MVL, V(1,p), 1,
Expand All @@ -1002,14 +1011,17 @@ SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
* .. have to use modified Gram-Schmidt like transformation
CALL ZCOPY( M, A( 1, p ), 1,
$ CWORK(N+1), 1 )
CALL ZLASCL( 'G', 0, 0, AAPP, ONE, M,
CALL ZLASCL( 'G', 0, 0, AAPP, ONE,
$ M,
$ 1, CWORK(N+1), LDA,
$ IERR )
CALL ZLASCL( 'G', 0, 0, AAQQ, ONE, M,
CALL ZLASCL( 'G', 0, 0, AAQQ, ONE,
$ M,
$ 1, A( 1, q ), LDA, IERR )
CALL ZAXPY( M, -AAPQ, CWORK(N+1), 1,
$ A( 1, q ), 1 )
CALL ZLASCL( 'G', 0, 0, ONE, AAQQ, M,
CALL ZLASCL( 'G', 0, 0, ONE, AAQQ,
$ M,
$ 1, A( 1, q ), LDA, IERR )
SVA( q ) = AAQQ*SQRT( MAX( ZERO,
$ ONE-AAPQ1*AAPQ1 ) )
Expand All @@ -1024,7 +1036,8 @@ SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
$ THEN
IF( ( AAQQ.LT.ROOTBIG ) .AND.
$ ( AAQQ.GT.ROOTSFMIN ) ) THEN
SVA( q ) = DZNRM2( M, A( 1, q ), 1 )
SVA( q ) = DZNRM2( M, A( 1, q ),
$ 1 )
ELSE
T = ZERO
AAQQ = ONE
Expand Down Expand Up @@ -1177,7 +1190,8 @@ SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
IF( ABS( THETA ).GT.BIGTHETA ) THEN
T = HALF / THETA
CS = ONE
CALL ZROT( M, A(1,p), 1, A(1,q), 1,
CALL ZROT( M, A(1,p), 1, A(1,q),
$ 1,
$ CS, CONJG(OMPQ)*T )
IF( RSVEC ) THEN
CALL ZROT( MVL, V(1,p), 1,
Expand All @@ -1204,7 +1218,8 @@ SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
AAPP = AAPP*SQRT( MAX( ZERO,
$ ONE-T*AQOAP*AAPQ1 ) )
*
CALL ZROT( M, A(1,p), 1, A(1,q), 1,
CALL ZROT( M, A(1,p), 1, A(1,q),
$ 1,
$ CS, CONJG(OMPQ)*SN )
IF( RSVEC ) THEN
CALL ZROT( MVL, V(1,p), 1,
Expand All @@ -1218,15 +1233,18 @@ SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
IF( AAPP.GT.AAQQ ) THEN
CALL ZCOPY( M, A( 1, p ), 1,
$ CWORK(N+1), 1 )
CALL ZLASCL( 'G', 0, 0, AAPP, ONE,
CALL ZLASCL( 'G', 0, 0, AAPP,
$ ONE,
$ M, 1, CWORK(N+1),LDA,
$ IERR )
CALL ZLASCL( 'G', 0, 0, AAQQ, ONE,
CALL ZLASCL( 'G', 0, 0, AAQQ,
$ ONE,
$ M, 1, A( 1, q ), LDA,
$ IERR )
CALL ZAXPY( M, -AAPQ, CWORK(N+1),
$ 1, A( 1, q ), 1 )
CALL ZLASCL( 'G', 0, 0, ONE, AAQQ,
CALL ZLASCL( 'G', 0, 0, ONE,
$ AAQQ,
$ M, 1, A( 1, q ), LDA,
$ IERR )
SVA( q ) = AAQQ*SQRT( MAX( ZERO,
Expand All @@ -1235,15 +1253,18 @@ SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
ELSE
CALL ZCOPY( M, A( 1, q ), 1,
$ CWORK(N+1), 1 )
CALL ZLASCL( 'G', 0, 0, AAQQ, ONE,
CALL ZLASCL( 'G', 0, 0, AAQQ,
$ ONE,
$ M, 1, CWORK(N+1),LDA,
$ IERR )
CALL ZLASCL( 'G', 0, 0, AAPP, ONE,
CALL ZLASCL( 'G', 0, 0, AAPP,
$ ONE,
$ M, 1, A( 1, p ), LDA,
$ IERR )
CALL ZAXPY( M, -CONJG(AAPQ),
$ CWORK(N+1), 1, A( 1, p ), 1 )
CALL ZLASCL( 'G', 0, 0, ONE, AAPP,
CALL ZLASCL( 'G', 0, 0, ONE,
$ AAPP,
$ M, 1, A( 1, p ), LDA,
$ IERR )
SVA( p ) = AAPP*SQRT( MAX( ZERO,
Expand All @@ -1259,7 +1280,8 @@ SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
$ THEN
IF( ( AAQQ.LT.ROOTBIG ) .AND.
$ ( AAQQ.GT.ROOTSFMIN ) ) THEN
SVA( q ) = DZNRM2( M, A( 1, q ), 1)
SVA( q ) = DZNRM2( M, A( 1, q ),
$ 1)
ELSE
T = ZERO
AAQQ = ONE
Expand Down Expand Up @@ -1401,7 +1423,8 @@ SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
IF( LSVEC .OR. UCTOL ) THEN
DO 1998 p = 1, N4
* CALL ZDSCAL( M, ONE / SVA( p ), A( 1, p ), 1 )
CALL ZLASCL( 'G',0,0, SVA(p), ONE, M, 1, A(1,p), M, IERR )
CALL ZLASCL( 'G',0,0, SVA(p), ONE, M, 1, A(1,p), M,
$ IERR )
1998 CONTINUE
END IF
*
Expand Down
24 changes: 14 additions & 10 deletions lapack-netlib/SRC/zhbevd.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,15 +5,13 @@
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download ZHBEVD + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhbevd.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhbevd.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhbevd.f">
*> [TXT]</a>
*> \endhtmlonly
*
* Definition:
* ===========
Expand Down Expand Up @@ -201,11 +199,13 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup complex16OTHEReigen
*> \ingroup hbevd
*
* =====================================================================
SUBROUTINE ZHBEVD( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
SUBROUTINE ZHBEVD( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ,
$ WORK,
$ LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO )
IMPLICIT NONE
*
* -- LAPACK driver routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
Expand Down Expand Up @@ -243,7 +243,8 @@ SUBROUTINE ZHBEVD( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
EXTERNAL LSAME, DLAMCH, ZLANHB
* ..
* .. External Subroutines ..
EXTERNAL DSCAL, DSTERF, XERBLA, ZGEMM, ZHBTRD, ZLACPY,
EXTERNAL DSCAL, DSTERF, XERBLA, ZGEMM, ZHBTRD,
$ ZLACPY,
$ ZLASCL, ZSTEDC
* ..
* .. Intrinsic Functions ..
Expand Down Expand Up @@ -289,7 +290,7 @@ SUBROUTINE ZHBEVD( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
*
IF( INFO.EQ.0 ) THEN
WORK( 1 ) = LWMIN
RWORK( 1 ) = LRWMIN
RWORK( 1 ) = DBLE( LRWMIN )
IWORK( 1 ) = LIWMIN
*
IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
Expand Down Expand Up @@ -342,9 +343,11 @@ SUBROUTINE ZHBEVD( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
END IF
IF( ISCALE.EQ.1 ) THEN
IF( LOWER ) THEN
CALL ZLASCL( 'B', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
CALL ZLASCL( 'B', KD, KD, ONE, SIGMA, N, N, AB, LDAB,
$ INFO )
ELSE
CALL ZLASCL( 'Q', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
CALL ZLASCL( 'Q', KD, KD, ONE, SIGMA, N, N, AB, LDAB,
$ INFO )
END IF
END IF
*
Expand All @@ -363,7 +366,8 @@ SUBROUTINE ZHBEVD( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
IF( .NOT.WANTZ ) THEN
CALL DSTERF( N, W, RWORK( INDE ), INFO )
ELSE
CALL ZSTEDC( 'I', N, W, RWORK( INDE ), WORK, N, WORK( INDWK2 ),
CALL ZSTEDC( 'I', N, W, RWORK( INDE ), WORK, N,
$ WORK( INDWK2 ),
$ LLWK2, RWORK( INDWRK ), LLRWK, IWORK, LIWORK,
$ INFO )
CALL ZGEMM( 'N', 'N', N, N, N, CONE, Z, LDZ, WORK, N, CZERO,
Expand All @@ -383,7 +387,7 @@ SUBROUTINE ZHBEVD( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
END IF
*
WORK( 1 ) = LWMIN
RWORK( 1 ) = LRWMIN
RWORK( 1 ) = DBLE( LRWMIN )
IWORK( 1 ) = LIWMIN
RETURN
*
Expand Down
Loading
Loading