Skip to content

Commit f1f36c0

Browse files
martin-frbgmartin-frbg
authored
Merge pull request #5729 from martin-frbg/lapack1195
Fix truncation of large workspace values in ZHE routines (Reference-LAPACK PR 1195)
2 parents 9816062 + 37e189c commit f1f36c0

10 files changed

Lines changed: 196 additions & 130 deletions

File tree

‎lapack-netlib/SRC/zgesvj.f‎

Lines changed: 53 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -5,15 +5,13 @@
55
* Online html documentation available at
66
* http://www.netlib.org/lapack/explore-html/
77
*
8-
*> \htmlonly
98
*> Download ZGESVJ + dependencies
109
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgesvj.f">
1110
*> [TGZ]</a>
1211
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgesvj.f">
1312
*> [ZIP]</a>
1413
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgesvj.f">
1514
*> [TXT]</a>
16-
*> \endhtmlonly
1715
*
1816
* Definition:
1917
* ===========
@@ -101,7 +99,7 @@
10199
*> \param[in] M
102100
*> \verbatim
103101
*> M is INTEGER
104-
*> The number of rows of the input matrix A. 1/DLAMCH('E') > M >= 0.
102+
*> The number of rows of the input matrix A. 1/DLAMCH('E') >= M >= 0.
105103
*> \endverbatim
106104
*>
107105
*> \param[in] N
@@ -217,7 +215,7 @@
217215
*> LWORK >= 1, if MIN(M,N) = 0, and LWORK >= M+N, otherwise.
218216
*>
219217
*> If on entry LWORK = -1, then a workspace query is assumed and
220-
*> no computation is done; CWORK(1) is set to the minial (and optimal)
218+
*> no computation is done; CWORK(1) is set to the minimal (and optimal)
221219
*> length of CWORK.
222220
*> \endverbatim
223221
*>
@@ -258,7 +256,7 @@
258256
*> LRWORK >= 1, if MIN(M,N) = 0, and LRWORK >= MAX(6,N), otherwise.
259257
*>
260258
*> If on entry LRWORK = -1, then a workspace query is assumed and
261-
*> no computation is done; RWORK(1) is set to the minial (and optimal)
259+
*> no computation is done; RWORK(1) is set to the minimal (and optimal)
262260
*> length of RWORK.
263261
*> \endverbatim
264262
*>
@@ -414,7 +412,8 @@ SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
414412
* from BLAS
415413
EXTERNAL ZCOPY, ZROT, ZDSCAL, ZSWAP, ZAXPY
416414
* from LAPACK
417-
EXTERNAL DLASCL, ZLASCL, ZLASET, ZLASSQ, XERBLA
415+
EXTERNAL DLASCL, ZLASCL, ZLASET, ZLASSQ,
416+
$ XERBLA
418417
EXTERNAL ZGSVJ0, ZGSVJ1
419418
* ..
420419
* .. Executable Statements ..
@@ -440,9 +439,13 @@ SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
440439
LQUERY = ( LWORK.EQ.-1 ) .OR. ( LRWORK.EQ.-1 )
441440
IF( .NOT.( UPPER .OR. LOWER .OR. LSAME( JOBA, 'G' ) ) ) THEN
442441
INFO = -1
443-
ELSE IF( .NOT.( LSVEC .OR. UCTOL .OR. LSAME( JOBU, 'N' ) ) ) THEN
442+
ELSE IF( .NOT.( LSVEC .OR.
443+
$ UCTOL .OR.
444+
$ LSAME( JOBU, 'N' ) ) ) THEN
444445
INFO = -2
445-
ELSE IF( .NOT.( RSVEC .OR. APPLV .OR. LSAME( JOBV, 'N' ) ) ) THEN
446+
ELSE IF( .NOT.( RSVEC .OR.
447+
$ APPLV .OR.
448+
$ LSAME( JOBV, 'N' ) ) ) THEN
446449
INFO = -3
447450
ELSE IF( M.LT.0 ) THEN
448451
INFO = -4
@@ -455,7 +458,7 @@ SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
455458
ELSE IF( ( RSVEC .AND. ( LDV.LT.N ) ) .OR.
456459
$ ( APPLV .AND. ( LDV.LT.MV ) ) ) THEN
457460
INFO = -11
458-
ELSE IF( UCTOL .AND. ( RWORK( 1 ).LE.ONE ) ) THEN
461+
ELSE IF( UCTOL .AND. ( RWORK( 1 ).LT.ONE ) ) THEN
459462
INFO = -12
460463
ELSE IF( LWORK.LT.LWMIN .AND. ( .NOT.LQUERY ) ) THEN
461464
INFO = -13
@@ -471,7 +474,7 @@ SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
471474
RETURN
472475
ELSE IF( LQUERY ) THEN
473476
CWORK( 1 ) = LWMIN
474-
RWORK( 1 ) = LRWMIN
477+
RWORK( 1 ) = DBLE( LRWMIN )
475478
RETURN
476479
END IF
477480
*
@@ -785,7 +788,8 @@ SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
785788
$ V( N4*q+1, N4+1 ), LDV, EPSLN, SFMIN, TOL, 1,
786789
$ CWORK( N+1 ), LWORK-N, IERR )
787790
*
788-
CALL ZGSVJ0( JOBV, M, N4, A, LDA, CWORK, SVA, MVL, V, LDV,
791+
CALL ZGSVJ0( JOBV, M, N4, A, LDA, CWORK, SVA, MVL, V,
792+
$ LDV,
789793
$ EPSLN, SFMIN, TOL, 1, CWORK( N+1 ), LWORK-N,
790794
$ IERR )
791795
*
@@ -797,16 +801,19 @@ SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
797801
ELSE IF( UPPER ) THEN
798802
*
799803
*
800-
CALL ZGSVJ0( JOBV, N4, N4, A, LDA, CWORK, SVA, MVL, V, LDV,
804+
CALL ZGSVJ0( JOBV, N4, N4, A, LDA, CWORK, SVA, MVL, V,
805+
$ LDV,
801806
$ EPSLN, SFMIN, TOL, 2, CWORK( N+1 ), LWORK-N,
802807
$ IERR )
803808
*
804-
CALL ZGSVJ0( JOBV, N2, N4, A( 1, N4+1 ), LDA, CWORK( N4+1 ),
809+
CALL ZGSVJ0( JOBV, N2, N4, A( 1, N4+1 ), LDA,
810+
$ CWORK( N4+1 ),
805811
$ SVA( N4+1 ), MVL, V( N4*q+1, N4+1 ), LDV,
806812
$ EPSLN, SFMIN, TOL, 1, CWORK( N+1 ), LWORK-N,
807813
$ IERR )
808814
*
809-
CALL ZGSVJ1( JOBV, N2, N2, N4, A, LDA, CWORK, SVA, MVL, V,
815+
CALL ZGSVJ1( JOBV, N2, N2, N4, A, LDA, CWORK, SVA, MVL,
816+
$ V,
810817
$ LDV, EPSLN, SFMIN, TOL, 1, CWORK( N+1 ),
811818
$ LWORK-N, IERR )
812819
*
@@ -960,7 +967,8 @@ SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
960967
T = HALF / THETA
961968
CS = ONE
962969

963-
CALL ZROT( M, A(1,p), 1, A(1,q), 1,
970+
CALL ZROT( M, A(1,p), 1, A(1,q),
971+
$ 1,
964972
$ CS, CONJG(OMPQ)*T )
965973
IF ( RSVEC ) THEN
966974
CALL ZROT( MVL, V(1,p), 1,
@@ -989,7 +997,8 @@ SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
989997
AAPP = AAPP*SQRT( MAX( ZERO,
990998
$ ONE-T*AQOAP*AAPQ1 ) )
991999
*
992-
CALL ZROT( M, A(1,p), 1, A(1,q), 1,
1000+
CALL ZROT( M, A(1,p), 1, A(1,q),
1001+
$ 1,
9931002
$ CS, CONJG(OMPQ)*SN )
9941003
IF ( RSVEC ) THEN
9951004
CALL ZROT( MVL, V(1,p), 1,
@@ -1002,14 +1011,17 @@ SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
10021011
* .. have to use modified Gram-Schmidt like transformation
10031012
CALL ZCOPY( M, A( 1, p ), 1,
10041013
$ CWORK(N+1), 1 )
1005-
CALL ZLASCL( 'G', 0, 0, AAPP, ONE, M,
1014+
CALL ZLASCL( 'G', 0, 0, AAPP, ONE,
1015+
$ M,
10061016
$ 1, CWORK(N+1), LDA,
10071017
$ IERR )
1008-
CALL ZLASCL( 'G', 0, 0, AAQQ, ONE, M,
1018+
CALL ZLASCL( 'G', 0, 0, AAQQ, ONE,
1019+
$ M,
10091020
$ 1, A( 1, q ), LDA, IERR )
10101021
CALL ZAXPY( M, -AAPQ, CWORK(N+1), 1,
10111022
$ A( 1, q ), 1 )
1012-
CALL ZLASCL( 'G', 0, 0, ONE, AAQQ, M,
1023+
CALL ZLASCL( 'G', 0, 0, ONE, AAQQ,
1024+
$ M,
10131025
$ 1, A( 1, q ), LDA, IERR )
10141026
SVA( q ) = AAQQ*SQRT( MAX( ZERO,
10151027
$ ONE-AAPQ1*AAPQ1 ) )
@@ -1024,7 +1036,8 @@ SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
10241036
$ THEN
10251037
IF( ( AAQQ.LT.ROOTBIG ) .AND.
10261038
$ ( AAQQ.GT.ROOTSFMIN ) ) THEN
1027-
SVA( q ) = DZNRM2( M, A( 1, q ), 1 )
1039+
SVA( q ) = DZNRM2( M, A( 1, q ),
1040+
$ 1 )
10281041
ELSE
10291042
T = ZERO
10301043
AAQQ = ONE
@@ -1177,7 +1190,8 @@ SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
11771190
IF( ABS( THETA ).GT.BIGTHETA ) THEN
11781191
T = HALF / THETA
11791192
CS = ONE
1180-
CALL ZROT( M, A(1,p), 1, A(1,q), 1,
1193+
CALL ZROT( M, A(1,p), 1, A(1,q),
1194+
$ 1,
11811195
$ CS, CONJG(OMPQ)*T )
11821196
IF( RSVEC ) THEN
11831197
CALL ZROT( MVL, V(1,p), 1,
@@ -1204,7 +1218,8 @@ SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
12041218
AAPP = AAPP*SQRT( MAX( ZERO,
12051219
$ ONE-T*AQOAP*AAPQ1 ) )
12061220
*
1207-
CALL ZROT( M, A(1,p), 1, A(1,q), 1,
1221+
CALL ZROT( M, A(1,p), 1, A(1,q),
1222+
$ 1,
12081223
$ CS, CONJG(OMPQ)*SN )
12091224
IF( RSVEC ) THEN
12101225
CALL ZROT( MVL, V(1,p), 1,
@@ -1218,15 +1233,18 @@ SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
12181233
IF( AAPP.GT.AAQQ ) THEN
12191234
CALL ZCOPY( M, A( 1, p ), 1,
12201235
$ CWORK(N+1), 1 )
1221-
CALL ZLASCL( 'G', 0, 0, AAPP, ONE,
1236+
CALL ZLASCL( 'G', 0, 0, AAPP,
1237+
$ ONE,
12221238
$ M, 1, CWORK(N+1),LDA,
12231239
$ IERR )
1224-
CALL ZLASCL( 'G', 0, 0, AAQQ, ONE,
1240+
CALL ZLASCL( 'G', 0, 0, AAQQ,
1241+
$ ONE,
12251242
$ M, 1, A( 1, q ), LDA,
12261243
$ IERR )
12271244
CALL ZAXPY( M, -AAPQ, CWORK(N+1),
12281245
$ 1, A( 1, q ), 1 )
1229-
CALL ZLASCL( 'G', 0, 0, ONE, AAQQ,
1246+
CALL ZLASCL( 'G', 0, 0, ONE,
1247+
$ AAQQ,
12301248
$ M, 1, A( 1, q ), LDA,
12311249
$ IERR )
12321250
SVA( q ) = AAQQ*SQRT( MAX( ZERO,
@@ -1235,15 +1253,18 @@ SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
12351253
ELSE
12361254
CALL ZCOPY( M, A( 1, q ), 1,
12371255
$ CWORK(N+1), 1 )
1238-
CALL ZLASCL( 'G', 0, 0, AAQQ, ONE,
1256+
CALL ZLASCL( 'G', 0, 0, AAQQ,
1257+
$ ONE,
12391258
$ M, 1, CWORK(N+1),LDA,
12401259
$ IERR )
1241-
CALL ZLASCL( 'G', 0, 0, AAPP, ONE,
1260+
CALL ZLASCL( 'G', 0, 0, AAPP,
1261+
$ ONE,
12421262
$ M, 1, A( 1, p ), LDA,
12431263
$ IERR )
12441264
CALL ZAXPY( M, -CONJG(AAPQ),
12451265
$ CWORK(N+1), 1, A( 1, p ), 1 )
1246-
CALL ZLASCL( 'G', 0, 0, ONE, AAPP,
1266+
CALL ZLASCL( 'G', 0, 0, ONE,
1267+
$ AAPP,
12471268
$ M, 1, A( 1, p ), LDA,
12481269
$ IERR )
12491270
SVA( p ) = AAPP*SQRT( MAX( ZERO,
@@ -1259,7 +1280,8 @@ SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
12591280
$ THEN
12601281
IF( ( AAQQ.LT.ROOTBIG ) .AND.
12611282
$ ( AAQQ.GT.ROOTSFMIN ) ) THEN
1262-
SVA( q ) = DZNRM2( M, A( 1, q ), 1)
1283+
SVA( q ) = DZNRM2( M, A( 1, q ),
1284+
$ 1)
12631285
ELSE
12641286
T = ZERO
12651287
AAQQ = ONE
@@ -1401,7 +1423,8 @@ SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
14011423
IF( LSVEC .OR. UCTOL ) THEN
14021424
DO 1998 p = 1, N4
14031425
* CALL ZDSCAL( M, ONE / SVA( p ), A( 1, p ), 1 )
1404-
CALL ZLASCL( 'G',0,0, SVA(p), ONE, M, 1, A(1,p), M, IERR )
1426+
CALL ZLASCL( 'G',0,0, SVA(p), ONE, M, 1, A(1,p), M,
1427+
$ IERR )
14051428
1998 CONTINUE
14061429
END IF
14071430
*

‎lapack-netlib/SRC/zhbevd.f‎

Lines changed: 14 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -5,15 +5,13 @@
55
* Online html documentation available at
66
* http://www.netlib.org/lapack/explore-html/
77
*
8-
*> \htmlonly
98
*> Download ZHBEVD + dependencies
109
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhbevd.f">
1110
*> [TGZ]</a>
1211
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhbevd.f">
1312
*> [ZIP]</a>
1413
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhbevd.f">
1514
*> [TXT]</a>
16-
*> \endhtmlonly
1715
*
1816
* Definition:
1917
* ===========
@@ -201,11 +199,13 @@
201199
*> \author Univ. of Colorado Denver
202200
*> \author NAG Ltd.
203201
*
204-
*> \ingroup complex16OTHEReigen
202+
*> \ingroup hbevd
205203
*
206204
* =====================================================================
207-
SUBROUTINE ZHBEVD( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
205+
SUBROUTINE ZHBEVD( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ,
206+
$ WORK,
208207
$ LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO )
208+
IMPLICIT NONE
209209
*
210210
* -- LAPACK driver routine --
211211
* -- LAPACK is a software package provided by Univ. of Tennessee, --
@@ -243,7 +243,8 @@ SUBROUTINE ZHBEVD( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
243243
EXTERNAL LSAME, DLAMCH, ZLANHB
244244
* ..
245245
* .. External Subroutines ..
246-
EXTERNAL DSCAL, DSTERF, XERBLA, ZGEMM, ZHBTRD, ZLACPY,
246+
EXTERNAL DSCAL, DSTERF, XERBLA, ZGEMM, ZHBTRD,
247+
$ ZLACPY,
247248
$ ZLASCL, ZSTEDC
248249
* ..
249250
* .. Intrinsic Functions ..
@@ -289,7 +290,7 @@ SUBROUTINE ZHBEVD( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
289290
*
290291
IF( INFO.EQ.0 ) THEN
291292
WORK( 1 ) = LWMIN
292-
RWORK( 1 ) = LRWMIN
293+
RWORK( 1 ) = DBLE( LRWMIN )
293294
IWORK( 1 ) = LIWMIN
294295
*
295296
IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
@@ -342,9 +343,11 @@ SUBROUTINE ZHBEVD( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
342343
END IF
343344
IF( ISCALE.EQ.1 ) THEN
344345
IF( LOWER ) THEN
345-
CALL ZLASCL( 'B', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
346+
CALL ZLASCL( 'B', KD, KD, ONE, SIGMA, N, N, AB, LDAB,
347+
$ INFO )
346348
ELSE
347-
CALL ZLASCL( 'Q', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
349+
CALL ZLASCL( 'Q', KD, KD, ONE, SIGMA, N, N, AB, LDAB,
350+
$ INFO )
348351
END IF
349352
END IF
350353
*
@@ -363,7 +366,8 @@ SUBROUTINE ZHBEVD( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
363366
IF( .NOT.WANTZ ) THEN
364367
CALL DSTERF( N, W, RWORK( INDE ), INFO )
365368
ELSE
366-
CALL ZSTEDC( 'I', N, W, RWORK( INDE ), WORK, N, WORK( INDWK2 ),
369+
CALL ZSTEDC( 'I', N, W, RWORK( INDE ), WORK, N,
370+
$ WORK( INDWK2 ),
367371
$ LLWK2, RWORK( INDWRK ), LLRWK, IWORK, LIWORK,
368372
$ INFO )
369373
CALL ZGEMM( 'N', 'N', N, N, N, CONE, Z, LDZ, WORK, N, CZERO,
@@ -383,7 +387,7 @@ SUBROUTINE ZHBEVD( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
383387
END IF
384388
*
385389
WORK( 1 ) = LWMIN
386-
RWORK( 1 ) = LRWMIN
390+
RWORK( 1 ) = DBLE( LRWMIN )
387391
IWORK( 1 ) = LIWMIN
388392
RETURN
389393
*

0 commit comments

Comments
 (0)