diff --git a/lapack-netlib/SRC/dgejsv.f b/lapack-netlib/SRC/dgejsv.f
index 1db85e9c26..bec88dfd15 100644
--- a/lapack-netlib/SRC/dgejsv.f
+++ b/lapack-netlib/SRC/dgejsv.f
@@ -5,7 +5,6 @@
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
-*> \htmlonly
*> Download DGEJSV + dependencies
*>
*> [TGZ]
@@ -13,7 +12,6 @@
*> [ZIP]
*>
*> [TXT]
-*> \endhtmlonly
*
* Definition:
* ===========
@@ -21,9 +19,9 @@
* SUBROUTINE DGEJSV( 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 ..
@@ -391,7 +389,7 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
-*> \ingroup doubleGEsing
+*> \ingroup gejsv
*
*> \par Further Details:
* =====================
@@ -473,13 +471,13 @@
SUBROUTINE DGEJSV( 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 ..
@@ -498,7 +496,7 @@ SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
* .. Local Scalars ..
DOUBLE PRECISION 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,
@@ -514,7 +512,8 @@ SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
EXTERNAL IDAMAX, LSAME, DLAMCH, DNRM2
* ..
* .. External Subroutines ..
- EXTERNAL DCOPY, DGELQF, DGEQP3, DGEQRF, DLACPY, DLASCL,
+ EXTERNAL DCOPY, DGELQF, DGEQP3, DGEQRF, DLACPY,
+ $ DLASCL,
$ DLASET, DLASSQ, DLASWP, DORGQR, DORMLQ,
$ DORMQR, DPOCON, DSCAL, DSWAP, DTRSM, XERBLA
*
@@ -629,7 +628,9 @@ SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
RETURN
END IF
AAQQ = DSQRT(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.
@@ -692,15 +693,20 @@ SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
CALL DLACPY( '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 DGEQRF( M, N, U,LDU, WORK, WORK(N+1),LWORK-N,IERR )
- CALL DORGQR( M,N1,1, U,LDU,WORK,WORK(N+1),LWORK-N,IERR )
+ CALL DGEQRF( M, N, U,LDU, WORK, WORK(N+1),LWORK-N,
+ $ IERR )
+ CALL DORGQR( M,N1,1, U,LDU,WORK,WORK(N+1),LWORK-N,
+ $ IERR )
CALL DCOPY( 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
@@ -1115,7 +1121,8 @@ SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
1949 CONTINUE
1947 CONTINUE
ELSE
- CALL DLASET( 'U', NR-1, NR-1, ZERO, ZERO, A(1,2), LDA )
+ CALL DLASET( 'U', NR-1, NR-1, ZERO, ZERO, A(1,2),
+ $ LDA )
END IF
*
* .. and one-sided Jacobi rotations are started on a lower
@@ -1139,7 +1146,8 @@ SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
DO 1998 p = 1, NR
CALL DCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 )
1998 CONTINUE
- CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
+ CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2),
+ $ LDV )
*
CALL DGESVJ( 'L','U','N', N, NR, V,LDV, SVA, NR, A,LDA,
$ WORK, LWORK, INFO )
@@ -1151,25 +1159,32 @@ SUBROUTINE DGEJSV( 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 DLASET( 'Lower', NR-1, NR-1, ZERO, ZERO, A(2,1), LDA )
- CALL DGELQF( NR, N, A, LDA, WORK, WORK(N+1), LWORK-N, IERR)
+ CALL DLASET( 'Lower', NR-1, NR-1, ZERO, ZERO, A(2,1),
+ $ LDA )
+ CALL DGELQF( NR, N, A, LDA, WORK, WORK(N+1), LWORK-N,
+ $ IERR)
CALL DLACPY( 'Lower', NR, NR, A, LDA, V, LDV )
- CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
+ CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2),
+ $ LDV )
CALL DGEQRF( NR, NR, V, LDV, WORK(N+1), WORK(2*N+1),
$ LWORK-2*N, IERR )
DO 8998 p = 1, NR
CALL DCOPY( NR-p+1, V(p,p), LDV, V(p,p), 1 )
8998 CONTINUE
- CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
+ CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2),
+ $ LDV )
*
CALL DGESVJ( 'Lower', 'U','N', NR, NR, V,LDV, SVA, NR, U,
- $ LDU, WORK(N+1), LWORK, INFO )
+ $ LDU, WORK(N+1), LWORK-N, INFO )
SCALEM = WORK(N+1)
NUMRANK = IDNINT(WORK(N+2))
IF ( NR .LT. N ) THEN
- CALL DLASET( 'A',N-NR, NR, ZERO,ZERO, V(NR+1,1), LDV )
- CALL DLASET( 'A',NR, N-NR, ZERO,ZERO, V(1,NR+1), LDV )
- CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE, V(NR+1,NR+1), LDV )
+ CALL DLASET( 'A',N-NR, NR, ZERO,ZERO, V(NR+1,1),
+ $ LDV )
+ CALL DLASET( 'A',NR, N-NR, ZERO,ZERO, V(1,NR+1),
+ $ LDV )
+ CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE, V(NR+1,NR+1),
+ $ LDV )
END IF
*
CALL DORMLQ( 'Left', 'Transpose', N, N, NR, A, LDA, WORK,
@@ -1213,8 +1228,10 @@ SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
IF ( NR .LT. M ) THEN
CALL DLASET( 'A', M-NR, NR,ZERO, ZERO, U(NR+1,1), LDU )
IF ( NR .LT. N1 ) THEN
- CALL DLASET( 'A',NR, N1-NR, ZERO, ZERO, U(1,NR+1), LDU )
- CALL DLASET( 'A',M-NR,N1-NR,ZERO,ONE,U(NR+1,NR+1), LDU )
+ CALL DLASET( 'A',NR, N1-NR, ZERO, ZERO, U(1,NR+1),
+ $ LDU )
+ CALL DLASET( 'A',M-NR,N1-NR,ZERO,ONE,U(NR+1,NR+1),
+ $ LDU )
END IF
END IF
*
@@ -1276,7 +1293,8 @@ SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
2968 CONTINUE
2969 CONTINUE
ELSE
- CALL DLASET( 'U', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
+ CALL DLASET( 'U', NR-1, NR-1, ZERO, ZERO, V(1,2),
+ $ LDV )
END IF
*
* Estimate the row scaled condition number of R1
@@ -1432,7 +1450,8 @@ SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
* equation is Q2*V2 = the product of the Jacobi rotations
* used in DGESVJ, premultiplied with the orthogonal matrix
* from the second QR factorization.
- CALL DTRSM( 'L','U','N','N', NR,NR,ONE, A,LDA, V,LDV )
+ CALL DTRSM( '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
@@ -1443,7 +1462,8 @@ SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
IF ( NR .LT. N ) THEN
CALL DLASET('A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV)
CALL DLASET('A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV)
- CALL DLASET('A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV)
+ CALL DLASET('A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),
+ $ LDV)
END IF
CALL DORMQR('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)
@@ -1457,7 +1477,8 @@ SUBROUTINE DGEJSV( 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 DGESVJ( 'L', 'U', 'N', NR, NR, V, LDV, SVA, NR, U,
+ CALL DGESVJ( '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 = IDNINT(WORK(2*N+N*NR+NR+2))
@@ -1465,7 +1486,8 @@ SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
CALL DCOPY( NR, V(1,p), 1, U(1,p), 1 )
CALL DSCAL( NR, SVA(p), U(1,p), 1 )
3870 CONTINUE
- CALL DTRSM('L','U','N','N',NR,NR,ONE,WORK(2*N+1),N,U,LDU)
+ CALL DTRSM('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
@@ -1478,7 +1500,8 @@ SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
IF ( NR .LT. N ) THEN
CALL DLASET( 'A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV )
CALL DLASET( 'A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV )
- CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV )
+ CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),
+ $ LDV )
END IF
CALL DORMQR( '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 )
@@ -1494,14 +1517,16 @@ SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
* defense ensures that DGEJSV completes the task.
* Compute the full SVD of L3 using DGESVJ with explicit
* accumulation of Jacobi rotations.
- CALL DGESVJ( 'L', 'U', 'V', NR, NR, V, LDV, SVA, NR, U,
+ CALL DGESVJ( '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 = IDNINT(WORK(2*N+N*NR+NR+2))
IF ( NR .LT. N ) THEN
CALL DLASET( 'A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV )
CALL DLASET( 'A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV )
- CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV )
+ CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),
+ $ LDV )
END IF
CALL DORMQR( '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 )
@@ -1539,10 +1564,12 @@ SUBROUTINE DGEJSV( 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 DLASET( 'A', M-NR, NR, ZERO, ZERO, U(NR+1,1), LDU )
+ CALL DLASET( 'A', M-NR, NR, ZERO, ZERO, U(NR+1,1),
+ $ LDU )
IF ( NR .LT. N1 ) THEN
CALL DLASET('A',NR,N1-NR,ZERO,ZERO,U(1,NR+1),LDU)
- CALL DLASET('A',M-NR,N1-NR,ZERO,ONE,U(NR+1,NR+1),LDU)
+ CALL DLASET('A',M-NR,N1-NR,ZERO,ONE,U(NR+1,NR+1),
+ $ LDU)
END IF
END IF
*
@@ -1611,8 +1638,10 @@ SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
IF ( N .LT. M ) THEN
CALL DLASET( 'A', M-N, N, ZERO, ZERO, U(N+1,1), LDU )
IF ( N .LT. N1 ) THEN
- CALL DLASET( 'A',N, N1-N, ZERO, ZERO, U(1,N+1),LDU )
- CALL DLASET( 'A',M-N,N1-N, ZERO, ONE,U(N+1,N+1),LDU )
+ CALL DLASET( 'A',N, N1-N, ZERO, ZERO, U(1,N+1),
+ $ LDU )
+ CALL DLASET( 'A',M-N,N1-N, ZERO, ONE,U(N+1,N+1),
+ $ LDU )
END IF
END IF
CALL DORMQR( 'Left', 'No Tr', M, N1, N, A, LDA, WORK, U,
@@ -1719,8 +1748,10 @@ SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
IF ( NR .LT. M ) THEN
CALL DLASET( 'A', M-NR, NR, ZERO, ZERO, U(NR+1,1), LDU )
IF ( NR .LT. N1 ) THEN
- CALL DLASET( 'A',NR, N1-NR, ZERO, ZERO, U(1,NR+1),LDU )
- CALL DLASET( 'A',M-NR,N1-NR, ZERO, ONE,U(NR+1,NR+1),LDU )
+ CALL DLASET( 'A',NR, N1-NR, ZERO, ZERO, U(1,NR+1),
+ $ LDU )
+ CALL DLASET( 'A',M-NR,N1-NR, ZERO, ONE,U(NR+1,NR+1),
+ $ LDU )
END IF
END IF
*
@@ -1745,7 +1776,8 @@ SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
* Undo scaling, if necessary (and possible)
*
IF ( USCAL2 .LE. (BIG/SVA(1))*USCAL1 ) THEN
- CALL DLASCL( 'G', 0, 0, USCAL1, USCAL2, NR, 1, SVA, N, IERR )
+ CALL DLASCL( 'G', 0, 0, USCAL1, USCAL2, NR, 1, SVA, N,
+ $ IERR )
USCAL1 = ONE
USCAL2 = ONE
END IF