Skip to content

Fix wrong JOBA type in SGESVJ call of SGEJSV (Reference-LAPACK PR 1269) #5848

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

Open
martin-frbg wants to merge 1 commit into
base: develop
Choose a base branch
from
+70 −38
Open

Fix wrong JOBA type in SGESVJ call of SGEJSV (Reference-LAPACK PR 1269) #5848

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
108 changes: 70 additions & 38 deletions lapack-netlib/SRC/sgejsv.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,25 +5,23 @@
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download SGEJSV + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgejsv.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgejsv.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgejsv.f">
*> [TXT]</a>
*> \endhtmlonly
*
* Definition:
* ===========
*
* SUBROUTINE SGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
* M, N, A, LDA, SVA, U, LDU, V, LDV,
* WORK, LWORK, IWORK, INFO )
* IMPLICIT NONE
*
* .. Scalar Arguments ..
* IMPLICIT NONE
* INTEGER INFO, LDA, LDU, LDV, LWORK, M, N
* ..
* .. Array Arguments ..
Expand Down Expand Up @@ -391,7 +389,7 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup realGEsing
*> \ingroup gejsv
*
*> \par Further Details:
* =====================
Expand Down Expand Up @@ -473,13 +471,13 @@
SUBROUTINE SGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
$ M, N, A, LDA, SVA, U, LDU, V, LDV,
$ WORK, LWORK, IWORK, INFO )
IMPLICIT NONE
*
* -- LAPACK computational routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
IMPLICIT NONE
INTEGER INFO, LDA, LDU, LDV, LWORK, M, N
* ..
* .. Array Arguments ..
Expand All @@ -498,7 +496,7 @@ SUBROUTINE SGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
* .. Local Scalars ..
REAL AAPP, AAQQ, AATMAX, AATMIN, BIG, BIG1, COND_OK,
$ CONDR1, CONDR2, ENTRA, ENTRAT, EPSLN, MAXPRJ, SCALEM,
$ SCONDA, SFMIN, SMALL, TEMP1, USCAL1, USCAL2, XSC
$ SCONDA, SFMIN, SMALL, TEMP1, USCAL1, USCAL2, XSC, TBIG
INTEGER IERR, N1, NR, NUMRANK, p, q, WARNING
LOGICAL ALMORT, DEFR, ERREST, GOSCAL, JRACC, KILL, LSVEC,
$ L2ABER, L2KILL, L2PERT, L2RANK, L2TRAN,
Expand All @@ -514,7 +512,8 @@ SUBROUTINE SGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
EXTERNAL ISAMAX, LSAME, SLAMCH, SNRM2
* ..
* .. External Subroutines ..
EXTERNAL SCOPY, SGELQF, SGEQP3, SGEQRF, SLACPY, SLASCL,
EXTERNAL SCOPY, SGELQF, SGEQP3, SGEQRF, SLACPY,
$ SLASCL,
$ SLASET, SLASSQ, SLASWP, SORGQR, SORMLQ,
$ SORMQR, SPOCON, SSCAL, SSWAP, STRSM, XERBLA
*
Expand Down Expand Up @@ -629,7 +628,9 @@ SUBROUTINE SGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
RETURN
END IF
AAQQ = SQRT(AAQQ)
IF ( ( AAPP .LT. (BIG / AAQQ) ) .AND. NOSCAL ) THEN
TBIG = BIG
IF( AAQQ.GT.ONE ) TBIG = BIG / AAQQ
IF ( ( AAPP .LT. TBIG ) .AND. NOSCAL ) THEN
SVA(p) = AAPP * AAQQ
ELSE
NOSCAL = .FALSE.
Expand Down Expand Up @@ -692,15 +693,20 @@ SUBROUTINE SGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
CALL SLACPY( 'A', M, 1, A, LDA, U, LDU )
* computing all M left singular vectors of the M x 1 matrix
IF ( N1 .NE. N ) THEN
CALL SGEQRF( M, N, U,LDU, WORK, WORK(N+1),LWORK-N,IERR )
CALL SORGQR( M,N1,1, U,LDU,WORK,WORK(N+1),LWORK-N,IERR )
CALL SGEQRF( M, N, U,LDU, WORK, WORK(N+1),LWORK-N,
$ IERR )
CALL SORGQR( M,N1,1, U,LDU,WORK,WORK(N+1),LWORK-N,
$ IERR )
CALL SCOPY( M, A(1,1), 1, U(1,1), 1 )
END IF
END IF
IF ( RSVEC ) THEN
V(1,1) = ONE
END IF
IF ( SVA(1) .LT. (BIG*SCALEM) ) THEN
TBIG = BIG
IF ( SCALEM .EQ. ZERO ) SCALEM = ONE
IF ( SCALEM .LT. ONE) TBIG = BIG*SCALEM
IF ( SVA(1) .LT. TBIG ) THEN
SVA(1) = SVA(1) / SCALEM
SCALEM = ONE
END IF
Expand Down Expand Up @@ -1115,7 +1121,8 @@ SUBROUTINE SGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
1949 CONTINUE
1947 CONTINUE
ELSE
CALL SLASET( 'U', NR-1, NR-1, ZERO, ZERO, A(1,2), LDA )
CALL SLASET( 'U', NR-1, NR-1, ZERO, ZERO, A(1,2),
$ LDA )
END IF
*
* .. and one-sided Jacobi rotations are started on a lower
Expand All @@ -1139,7 +1146,8 @@ SUBROUTINE SGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
DO 1998 p = 1, NR
CALL SCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 )
1998 CONTINUE
CALL SLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
CALL SLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2),
$ LDV )
*
CALL SGESVJ( 'L','U','N', N, NR, V,LDV, SVA, NR, A,LDA,
$ WORK, LWORK, INFO )
Expand All @@ -1151,25 +1159,32 @@ SUBROUTINE SGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
* .. two more QR factorizations ( one QRF is not enough, two require
* accumulated product of Jacobi rotations, three are perfect )
*
CALL SLASET( 'Lower', NR-1, NR-1, ZERO, ZERO, A(2,1), LDA )
CALL SGELQF( NR, N, A, LDA, WORK, WORK(N+1), LWORK-N, IERR)
CALL SLASET( 'Lower', NR-1, NR-1, ZERO, ZERO, A(2,1),
$ LDA )
CALL SGELQF( NR, N, A, LDA, WORK, WORK(N+1), LWORK-N,
$ IERR)
CALL SLACPY( 'Lower', NR, NR, A, LDA, V, LDV )
CALL SLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
CALL SLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2),
$ LDV )
CALL SGEQRF( NR, NR, V, LDV, WORK(N+1), WORK(2*N+1),
$ LWORK-2*N, IERR )
DO 8998 p = 1, NR
CALL SCOPY( NR-p+1, V(p,p), LDV, V(p,p), 1 )
8998 CONTINUE
CALL SLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
CALL SLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2),
$ LDV )
*
CALL SGESVJ( 'Lower', 'U','N', NR, NR, V,LDV, SVA, NR, U,
$ LDU, WORK(N+1), LWORK-N, INFO )
SCALEM = WORK(N+1)
NUMRANK = NINT(WORK(N+2))
IF ( NR .LT. N ) THEN
CALL SLASET( 'A',N-NR, NR, ZERO,ZERO, V(NR+1,1), LDV )
CALL SLASET( 'A',NR, N-NR, ZERO,ZERO, V(1,NR+1), LDV )
CALL SLASET( 'A',N-NR,N-NR,ZERO,ONE, V(NR+1,NR+1), LDV )
CALL SLASET( 'A',N-NR, NR, ZERO,ZERO, V(NR+1,1),
$ LDV )
CALL SLASET( 'A',NR, N-NR, ZERO,ZERO, V(1,NR+1),
$ LDV )
CALL SLASET( 'A',N-NR,N-NR,ZERO,ONE, V(NR+1,NR+1),
$ LDV )
END IF
*
CALL SORMLQ( 'Left', 'Transpose', N, N, NR, A, LDA, WORK,
Expand Down Expand Up @@ -1213,8 +1228,10 @@ SUBROUTINE SGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
IF ( NR .LT. M ) THEN
CALL SLASET( 'A', M-NR, NR,ZERO, ZERO, U(NR+1,1), LDU )
IF ( NR .LT. N1 ) THEN
CALL SLASET( 'A',NR, N1-NR, ZERO, ZERO, U(1,NR+1), LDU )
CALL SLASET( 'A',M-NR,N1-NR,ZERO,ONE,U(NR+1,NR+1), LDU )
CALL SLASET( 'A',NR, N1-NR, ZERO, ZERO, U(1,NR+1),
$ LDU )
CALL SLASET( 'A',M-NR,N1-NR,ZERO,ONE,U(NR+1,NR+1),
$ LDU )
END IF
END IF
*
Expand Down Expand Up @@ -1276,7 +1293,8 @@ SUBROUTINE SGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
2968 CONTINUE
2969 CONTINUE
ELSE
CALL SLASET( 'U', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
CALL SLASET( 'U', NR-1, NR-1, ZERO, ZERO, V(1,2),
$ LDV )
END IF
*
* Estimate the row scaled condition number of R1
Expand Down Expand Up @@ -1432,7 +1450,8 @@ SUBROUTINE SGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
* equation is Q2*V2 = the product of the Jacobi rotations
* used in SGESVJ, premultiplied with the orthogonal matrix
* from the second QR factorization.
CALL STRSM( 'L','U','N','N', NR,NR,ONE, A,LDA, V,LDV )
CALL STRSM( 'L','U','N','N', NR,NR,ONE, A,LDA, V,
$ LDV )
ELSE
* .. R1 is well conditioned, but non-square. Transpose(R2)
* is inverted to get the product of the Jacobi rotations
Expand All @@ -1443,7 +1462,8 @@ SUBROUTINE SGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
IF ( NR .LT. N ) THEN
CALL SLASET('A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV)
CALL SLASET('A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV)
CALL SLASET('A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV)
CALL SLASET('A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),
$ LDV)
END IF
CALL SORMQR('L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1),
$ V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR)
Expand All @@ -1457,15 +1477,17 @@ SUBROUTINE SGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
* is Q3^T*V3 = the product of the Jacobi rotations (applied to
* the lower triangular L3 from the LQ factorization of
* R2=L3*Q3), pre-multiplied with the transposed Q3.
CALL SGESVJ( 'L', 'U', 'N', NR, NR, V, LDV, SVA, NR, U,
CALL SGESVJ( 'L', 'U', 'N', NR, NR, V, LDV, SVA, NR,
$ U,
$ LDU, WORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, INFO )
SCALEM = WORK(2*N+N*NR+NR+1)
NUMRANK = NINT(WORK(2*N+N*NR+NR+2))
DO 3870 p = 1, NR
CALL SCOPY( NR, V(1,p), 1, U(1,p), 1 )
CALL SSCAL( NR, SVA(p), U(1,p), 1 )
3870 CONTINUE
CALL STRSM('L','U','N','N',NR,NR,ONE,WORK(2*N+1),N,U,LDU)
CALL STRSM('L','U','N','N',NR,NR,ONE,WORK(2*N+1),N,U,
$ LDU)
* .. apply the permutation from the second QR factorization
DO 873 q = 1, NR
DO 872 p = 1, NR
Expand All @@ -1478,7 +1500,8 @@ SUBROUTINE SGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
IF ( NR .LT. N ) THEN
CALL SLASET( 'A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV )
CALL SLASET( 'A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV )
CALL SLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV )
CALL SLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),
$ LDV )
END IF
CALL SORMQR( 'L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1),
$ V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR )
Expand All @@ -1494,14 +1517,16 @@ SUBROUTINE SGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
* defense ensures that SGEJSV completes the task.
* Compute the full SVD of L3 using SGESVJ with explicit
* accumulation of Jacobi rotations.
CALL SGESVJ( 'L', 'U', 'V', NR, NR, V, LDV, SVA, NR, U,
CALL SGESVJ( 'L', 'U', 'V', NR, NR, V, LDV, SVA, NR,
$ U,
$ LDU, WORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, INFO )
SCALEM = WORK(2*N+N*NR+NR+1)
NUMRANK = NINT(WORK(2*N+N*NR+NR+2))
IF ( NR .LT. N ) THEN
CALL SLASET( 'A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV )
CALL SLASET( 'A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV )
CALL SLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV )
CALL SLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),
$ LDV )
END IF
CALL SORMQR( 'L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1),
$ V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR )
Expand Down Expand Up @@ -1539,10 +1564,12 @@ SUBROUTINE SGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
* At this moment, V contains the right singular vectors of A.
* Next, assemble the left singular vector matrix U (M x N).
IF ( NR .LT. M ) THEN
CALL SLASET( 'A', M-NR, NR, ZERO, ZERO, U(NR+1,1), LDU )
CALL SLASET( 'A', M-NR, NR, ZERO, ZERO, U(NR+1,1),
$ LDU )
IF ( NR .LT. N1 ) THEN
CALL SLASET('A',NR,N1-NR,ZERO,ZERO,U(1,NR+1),LDU)
CALL SLASET('A',M-NR,N1-NR,ZERO,ONE,U(NR+1,NR+1),LDU)
CALL SLASET('A',M-NR,N1-NR,ZERO,ONE,U(NR+1,NR+1),
$ LDU)
END IF
END IF
*
Expand Down Expand Up @@ -1611,8 +1638,10 @@ SUBROUTINE SGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
IF ( N .LT. M ) THEN
CALL SLASET( 'A', M-N, N, ZERO, ZERO, U(N+1,1), LDU )
IF ( N .LT. N1 ) THEN
CALL SLASET( 'A',N, N1-N, ZERO, ZERO, U(1,N+1),LDU )
CALL SLASET( 'A',M-N,N1-N, ZERO, ONE,U(N+1,N+1),LDU )
CALL SLASET( 'A',N, N1-N, ZERO, ZERO, U(1,N+1),
$ LDU )
CALL SLASET( 'A',M-N,N1-N, ZERO, ONE,U(N+1,N+1),
$ LDU )
END IF
END IF
CALL SORMQR( 'Left', 'No Tr', M, N1, N, A, LDA, WORK, U,
Expand Down Expand Up @@ -1682,7 +1711,7 @@ SUBROUTINE SGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
CALL SLASET('U', NR-1, NR-1, ZERO, ZERO, U(1,2), LDU )
END IF

CALL SGESVJ( 'L', 'U', 'V', NR, NR, U, LDU, SVA,
CALL SGESVJ( 'G', 'U', 'V', NR, NR, U, LDU, SVA,
$ N, V, LDV, WORK(2*N+N*NR+1), LWORK-2*N-N*NR, INFO )
SCALEM = WORK(2*N+N*NR+1)
NUMRANK = NINT(WORK(2*N+N*NR+2))
Expand Down Expand Up @@ -1719,8 +1748,10 @@ SUBROUTINE SGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
IF ( NR .LT. M ) THEN
CALL SLASET( 'A', M-NR, NR, ZERO, ZERO, U(NR+1,1), LDU )
IF ( NR .LT. N1 ) THEN
CALL SLASET( 'A',NR, N1-NR, ZERO, ZERO, U(1,NR+1),LDU )
CALL SLASET( 'A',M-NR,N1-NR, ZERO, ONE,U(NR+1,NR+1),LDU )
CALL SLASET( 'A',NR, N1-NR, ZERO, ZERO, U(1,NR+1),
$ LDU )
CALL SLASET( 'A',M-NR,N1-NR, ZERO, ONE,U(NR+1,NR+1),
$ LDU )
END IF
END IF
*
Expand All @@ -1745,7 +1776,8 @@ SUBROUTINE SGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
* Undo scaling, if necessary (and possible)
*
IF ( USCAL2 .LE. (BIG/SVA(1))*USCAL1 ) THEN
CALL SLASCL( 'G', 0, 0, USCAL1, USCAL2, NR, 1, SVA, N, IERR )
CALL SLASCL( 'G', 0, 0, USCAL1, USCAL2, NR, 1, SVA, N,
$ IERR )
USCAL1 = ONE
USCAL2 = ONE
END IF
Expand Down
Loading