Skip to content

Commit 1afb4c9

Browse files
martin-frbgmartin-frbg
authored
Merge pull request #5849 from martin-frbg/lapack1270
Fix wrong workspace in DGEJSV potentially corrupting memory in DGESVJ (Reference-LAPACK PR 1270)
2 parents 6d648c6 + 9f6c07c commit 1afb4c9

1 file changed

Lines changed: 70 additions & 38 deletions

File tree

‎lapack-netlib/SRC/dgejsv.f‎

Lines changed: 70 additions & 38 deletions
Original file line numberDiff line numberDiff line change
@@ -5,25 +5,23 @@
55
* Online html documentation available at
66
* http://www.netlib.org/lapack/explore-html/
77
*
8-
*> \htmlonly
98
*> Download DGEJSV + dependencies
109
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgejsv.f">
1110
*> [TGZ]</a>
1211
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgejsv.f">
1312
*> [ZIP]</a>
1413
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgejsv.f">
1514
*> [TXT]</a>
16-
*> \endhtmlonly
1715
*
1816
* Definition:
1917
* ===========
2018
*
2119
* SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
2220
* M, N, A, LDA, SVA, U, LDU, V, LDV,
2321
* WORK, LWORK, IWORK, INFO )
22+
* IMPLICIT NONE
2423
*
2524
* .. Scalar Arguments ..
26-
* IMPLICIT NONE
2725
* INTEGER INFO, LDA, LDU, LDV, LWORK, M, N
2826
* ..
2927
* .. Array Arguments ..
@@ -391,7 +389,7 @@
391389
*> \author Univ. of Colorado Denver
392390
*> \author NAG Ltd.
393391
*
394-
*> \ingroup doubleGEsing
392+
*> \ingroup gejsv
395393
*
396394
*> \par Further Details:
397395
* =====================
@@ -473,13 +471,13 @@
473471
SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
474472
$ M, N, A, LDA, SVA, U, LDU, V, LDV,
475473
$ WORK, LWORK, IWORK, INFO )
474+
IMPLICIT NONE
476475
*
477476
* -- LAPACK computational routine --
478477
* -- LAPACK is a software package provided by Univ. of Tennessee, --
479478
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
480479
*
481480
* .. Scalar Arguments ..
482-
IMPLICIT NONE
483481
INTEGER INFO, LDA, LDU, LDV, LWORK, M, N
484482
* ..
485483
* .. Array Arguments ..
@@ -498,7 +496,7 @@ SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
498496
* .. Local Scalars ..
499497
DOUBLE PRECISION AAPP, AAQQ, AATMAX, AATMIN, BIG, BIG1, COND_OK,
500498
$ CONDR1, CONDR2, ENTRA, ENTRAT, EPSLN, MAXPRJ, SCALEM,
501-
$ SCONDA, SFMIN, SMALL, TEMP1, USCAL1, USCAL2, XSC
499+
$ SCONDA, SFMIN, SMALL, TEMP1, USCAL1, USCAL2, XSC, TBIG
502500
INTEGER IERR, N1, NR, NUMRANK, p, q, WARNING
503501
LOGICAL ALMORT, DEFR, ERREST, GOSCAL, JRACC, KILL, LSVEC,
504502
$ L2ABER, L2KILL, L2PERT, L2RANK, L2TRAN,
@@ -514,7 +512,8 @@ SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
514512
EXTERNAL IDAMAX, LSAME, DLAMCH, DNRM2
515513
* ..
516514
* .. External Subroutines ..
517-
EXTERNAL DCOPY, DGELQF, DGEQP3, DGEQRF, DLACPY, DLASCL,
515+
EXTERNAL DCOPY, DGELQF, DGEQP3, DGEQRF, DLACPY,
516+
$ DLASCL,
518517
$ DLASET, DLASSQ, DLASWP, DORGQR, DORMLQ,
519518
$ DORMQR, DPOCON, DSCAL, DSWAP, DTRSM, XERBLA
520519
*
@@ -629,7 +628,9 @@ SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
629628
RETURN
630629
END IF
631630
AAQQ = DSQRT(AAQQ)
632-
IF ( ( AAPP .LT. (BIG / AAQQ) ) .AND. NOSCAL ) THEN
631+
TBIG = BIG
632+
IF( AAQQ.GT.ONE ) TBIG = BIG / AAQQ
633+
IF ( ( AAPP .LT. TBIG ) .AND. NOSCAL ) THEN
633634
SVA(p) = AAPP * AAQQ
634635
ELSE
635636
NOSCAL = .FALSE.
@@ -692,15 +693,20 @@ SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
692693
CALL DLACPY( 'A', M, 1, A, LDA, U, LDU )
693694
* computing all M left singular vectors of the M x 1 matrix
694695
IF ( N1 .NE. N ) THEN
695-
CALL DGEQRF( M, N, U,LDU, WORK, WORK(N+1),LWORK-N,IERR )
696-
CALL DORGQR( M,N1,1, U,LDU,WORK,WORK(N+1),LWORK-N,IERR )
696+
CALL DGEQRF( M, N, U,LDU, WORK, WORK(N+1),LWORK-N,
697+
$ IERR )
698+
CALL DORGQR( M,N1,1, U,LDU,WORK,WORK(N+1),LWORK-N,
699+
$ IERR )
697700
CALL DCOPY( M, A(1,1), 1, U(1,1), 1 )
698701
END IF
699702
END IF
700703
IF ( RSVEC ) THEN
701704
V(1,1) = ONE
702705
END IF
703-
IF ( SVA(1) .LT. (BIG*SCALEM) ) THEN
706+
TBIG = BIG
707+
IF ( SCALEM .EQ. ZERO ) SCALEM = ONE
708+
IF ( SCALEM .LT. ONE) TBIG = BIG*SCALEM
709+
IF ( SVA(1) .LT. TBIG ) THEN
704710
SVA(1) = SVA(1) / SCALEM
705711
SCALEM = ONE
706712
END IF
@@ -1115,7 +1121,8 @@ SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
11151121
1949 CONTINUE
11161122
1947 CONTINUE
11171123
ELSE
1118-
CALL DLASET( 'U', NR-1, NR-1, ZERO, ZERO, A(1,2), LDA )
1124+
CALL DLASET( 'U', NR-1, NR-1, ZERO, ZERO, A(1,2),
1125+
$ LDA )
11191126
END IF
11201127
*
11211128
* .. and one-sided Jacobi rotations are started on a lower
@@ -1139,7 +1146,8 @@ SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
11391146
DO 1998 p = 1, NR
11401147
CALL DCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 )
11411148
1998 CONTINUE
1142-
CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
1149+
CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2),
1150+
$ LDV )
11431151
*
11441152
CALL DGESVJ( 'L','U','N', N, NR, V,LDV, SVA, NR, A,LDA,
11451153
$ WORK, LWORK, INFO )
@@ -1151,25 +1159,32 @@ SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
11511159
* .. two more QR factorizations ( one QRF is not enough, two require
11521160
* accumulated product of Jacobi rotations, three are perfect )
11531161
*
1154-
CALL DLASET( 'Lower', NR-1, NR-1, ZERO, ZERO, A(2,1), LDA )
1155-
CALL DGELQF( NR, N, A, LDA, WORK, WORK(N+1), LWORK-N, IERR)
1162+
CALL DLASET( 'Lower', NR-1, NR-1, ZERO, ZERO, A(2,1),
1163+
$ LDA )
1164+
CALL DGELQF( NR, N, A, LDA, WORK, WORK(N+1), LWORK-N,
1165+
$ IERR)
11561166
CALL DLACPY( 'Lower', NR, NR, A, LDA, V, LDV )
1157-
CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
1167+
CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2),
1168+
$ LDV )
11581169
CALL DGEQRF( NR, NR, V, LDV, WORK(N+1), WORK(2*N+1),
11591170
$ LWORK-2*N, IERR )
11601171
DO 8998 p = 1, NR
11611172
CALL DCOPY( NR-p+1, V(p,p), LDV, V(p,p), 1 )
11621173
8998 CONTINUE
1163-
CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
1174+
CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2),
1175+
$ LDV )
11641176
*
11651177
CALL DGESVJ( 'Lower', 'U','N', NR, NR, V,LDV, SVA, NR, U,
1166-
$ LDU, WORK(N+1), LWORK, INFO )
1178+
$ LDU, WORK(N+1), LWORK-N, INFO )
11671179
SCALEM = WORK(N+1)
11681180
NUMRANK = IDNINT(WORK(N+2))
11691181
IF ( NR .LT. N ) THEN
1170-
CALL DLASET( 'A',N-NR, NR, ZERO,ZERO, V(NR+1,1), LDV )
1171-
CALL DLASET( 'A',NR, N-NR, ZERO,ZERO, V(1,NR+1), LDV )
1172-
CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE, V(NR+1,NR+1), LDV )
1182+
CALL DLASET( 'A',N-NR, NR, ZERO,ZERO, V(NR+1,1),
1183+
$ LDV )
1184+
CALL DLASET( 'A',NR, N-NR, ZERO,ZERO, V(1,NR+1),
1185+
$ LDV )
1186+
CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE, V(NR+1,NR+1),
1187+
$ LDV )
11731188
END IF
11741189
*
11751190
CALL DORMLQ( 'Left', 'Transpose', N, N, NR, A, LDA, WORK,
@@ -1213,8 +1228,10 @@ SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
12131228
IF ( NR .LT. M ) THEN
12141229
CALL DLASET( 'A', M-NR, NR,ZERO, ZERO, U(NR+1,1), LDU )
12151230
IF ( NR .LT. N1 ) THEN
1216-
CALL DLASET( 'A',NR, N1-NR, ZERO, ZERO, U(1,NR+1), LDU )
1217-
CALL DLASET( 'A',M-NR,N1-NR,ZERO,ONE,U(NR+1,NR+1), LDU )
1231+
CALL DLASET( 'A',NR, N1-NR, ZERO, ZERO, U(1,NR+1),
1232+
$ LDU )
1233+
CALL DLASET( 'A',M-NR,N1-NR,ZERO,ONE,U(NR+1,NR+1),
1234+
$ LDU )
12181235
END IF
12191236
END IF
12201237
*
@@ -1276,7 +1293,8 @@ SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
12761293
2968 CONTINUE
12771294
2969 CONTINUE
12781295
ELSE
1279-
CALL DLASET( 'U', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
1296+
CALL DLASET( 'U', NR-1, NR-1, ZERO, ZERO, V(1,2),
1297+
$ LDV )
12801298
END IF
12811299
*
12821300
* Estimate the row scaled condition number of R1
@@ -1432,7 +1450,8 @@ SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
14321450
* equation is Q2*V2 = the product of the Jacobi rotations
14331451
* used in DGESVJ, premultiplied with the orthogonal matrix
14341452
* from the second QR factorization.
1435-
CALL DTRSM( 'L','U','N','N', NR,NR,ONE, A,LDA, V,LDV )
1453+
CALL DTRSM( 'L','U','N','N', NR,NR,ONE, A,LDA, V,
1454+
$ LDV )
14361455
ELSE
14371456
* .. R1 is well conditioned, but non-square. Transpose(R2)
14381457
* is inverted to get the product of the Jacobi rotations
@@ -1443,7 +1462,8 @@ SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
14431462
IF ( NR .LT. N ) THEN
14441463
CALL DLASET('A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV)
14451464
CALL DLASET('A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV)
1446-
CALL DLASET('A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV)
1465+
CALL DLASET('A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),
1466+
$ LDV)
14471467
END IF
14481468
CALL DORMQR('L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1),
14491469
$ V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR)
@@ -1457,15 +1477,17 @@ SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
14571477
* is Q3^T*V3 = the product of the Jacobi rotations (applied to
14581478
* the lower triangular L3 from the LQ factorization of
14591479
* R2=L3*Q3), pre-multiplied with the transposed Q3.
1460-
CALL DGESVJ( 'L', 'U', 'N', NR, NR, V, LDV, SVA, NR, U,
1480+
CALL DGESVJ( 'L', 'U', 'N', NR, NR, V, LDV, SVA, NR,
1481+
$ U,
14611482
$ LDU, WORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, INFO )
14621483
SCALEM = WORK(2*N+N*NR+NR+1)
14631484
NUMRANK = IDNINT(WORK(2*N+N*NR+NR+2))
14641485
DO 3870 p = 1, NR
14651486
CALL DCOPY( NR, V(1,p), 1, U(1,p), 1 )
14661487
CALL DSCAL( NR, SVA(p), U(1,p), 1 )
14671488
3870 CONTINUE
1468-
CALL DTRSM('L','U','N','N',NR,NR,ONE,WORK(2*N+1),N,U,LDU)
1489+
CALL DTRSM('L','U','N','N',NR,NR,ONE,WORK(2*N+1),N,U,
1490+
$ LDU)
14691491
* .. apply the permutation from the second QR factorization
14701492
DO 873 q = 1, NR
14711493
DO 872 p = 1, NR
@@ -1478,7 +1500,8 @@ SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
14781500
IF ( NR .LT. N ) THEN
14791501
CALL DLASET( 'A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV )
14801502
CALL DLASET( 'A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV )
1481-
CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV )
1503+
CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),
1504+
$ LDV )
14821505
END IF
14831506
CALL DORMQR( 'L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1),
14841507
$ V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR )
@@ -1494,14 +1517,16 @@ SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
14941517
* defense ensures that DGEJSV completes the task.
14951518
* Compute the full SVD of L3 using DGESVJ with explicit
14961519
* accumulation of Jacobi rotations.
1497-
CALL DGESVJ( 'L', 'U', 'V', NR, NR, V, LDV, SVA, NR, U,
1520+
CALL DGESVJ( 'L', 'U', 'V', NR, NR, V, LDV, SVA, NR,
1521+
$ U,
14981522
$ LDU, WORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, INFO )
14991523
SCALEM = WORK(2*N+N*NR+NR+1)
15001524
NUMRANK = IDNINT(WORK(2*N+N*NR+NR+2))
15011525
IF ( NR .LT. N ) THEN
15021526
CALL DLASET( 'A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV )
15031527
CALL DLASET( 'A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV )
1504-
CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV )
1528+
CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),
1529+
$ LDV )
15051530
END IF
15061531
CALL DORMQR( 'L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1),
15071532
$ V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR )
@@ -1539,10 +1564,12 @@ SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
15391564
* At this moment, V contains the right singular vectors of A.
15401565
* Next, assemble the left singular vector matrix U (M x N).
15411566
IF ( NR .LT. M ) THEN
1542-
CALL DLASET( 'A', M-NR, NR, ZERO, ZERO, U(NR+1,1), LDU )
1567+
CALL DLASET( 'A', M-NR, NR, ZERO, ZERO, U(NR+1,1),
1568+
$ LDU )
15431569
IF ( NR .LT. N1 ) THEN
15441570
CALL DLASET('A',NR,N1-NR,ZERO,ZERO,U(1,NR+1),LDU)
1545-
CALL DLASET('A',M-NR,N1-NR,ZERO,ONE,U(NR+1,NR+1),LDU)
1571+
CALL DLASET('A',M-NR,N1-NR,ZERO,ONE,U(NR+1,NR+1),
1572+
$ LDU)
15461573
END IF
15471574
END IF
15481575
*
@@ -1611,8 +1638,10 @@ SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
16111638
IF ( N .LT. M ) THEN
16121639
CALL DLASET( 'A', M-N, N, ZERO, ZERO, U(N+1,1), LDU )
16131640
IF ( N .LT. N1 ) THEN
1614-
CALL DLASET( 'A',N, N1-N, ZERO, ZERO, U(1,N+1),LDU )
1615-
CALL DLASET( 'A',M-N,N1-N, ZERO, ONE,U(N+1,N+1),LDU )
1641+
CALL DLASET( 'A',N, N1-N, ZERO, ZERO, U(1,N+1),
1642+
$ LDU )
1643+
CALL DLASET( 'A',M-N,N1-N, ZERO, ONE,U(N+1,N+1),
1644+
$ LDU )
16161645
END IF
16171646
END IF
16181647
CALL DORMQR( 'Left', 'No Tr', M, N1, N, A, LDA, WORK, U,
@@ -1719,8 +1748,10 @@ SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
17191748
IF ( NR .LT. M ) THEN
17201749
CALL DLASET( 'A', M-NR, NR, ZERO, ZERO, U(NR+1,1), LDU )
17211750
IF ( NR .LT. N1 ) THEN
1722-
CALL DLASET( 'A',NR, N1-NR, ZERO, ZERO, U(1,NR+1),LDU )
1723-
CALL DLASET( 'A',M-NR,N1-NR, ZERO, ONE,U(NR+1,NR+1),LDU )
1751+
CALL DLASET( 'A',NR, N1-NR, ZERO, ZERO, U(1,NR+1),
1752+
$ LDU )
1753+
CALL DLASET( 'A',M-NR,N1-NR, ZERO, ONE,U(NR+1,NR+1),
1754+
$ LDU )
17241755
END IF
17251756
END IF
17261757
*
@@ -1745,7 +1776,8 @@ SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
17451776
* Undo scaling, if necessary (and possible)
17461777
*
17471778
IF ( USCAL2 .LE. (BIG/SVA(1))*USCAL1 ) THEN
1748-
CALL DLASCL( 'G', 0, 0, USCAL1, USCAL2, NR, 1, SVA, N, IERR )
1779+
CALL DLASCL( 'G', 0, 0, USCAL1, USCAL2, NR, 1, SVA, N,
1780+
$ IERR )
17491781
USCAL1 = ONE
17501782
USCAL2 = ONE
17511783
END IF

0 commit comments

Comments
 (0)