GOSCAL = .TRUE.
DO 1874 p = 1, N
AAPP = ZERO
- AAQQ = ZERO
+ AAQQ = ONE
CALL DLASSQ( M, A(1,p), 1, AAPP, AAQQ )
IF ( AAPP .GT. BIG ) THEN
INFO = - 9
IF ( L2TRAN ) THEN
DO 1950 p = 1, M
XSC = ZERO
- TEMP1 = ZERO
+ TEMP1 = ONE
CALL DLASSQ( N, A(p,1), LDA, XSC, TEMP1 )
* DLASSQ gets both the ell_2 and the ell_infinity norm
* in one pass through the vector
IF ( L2TRAN ) THEN
*
XSC = ZERO
- TEMP1 = ZERO
+ TEMP1 = ONE
CALL DLASSQ( N, SVA, 1, XSC, TEMP1 )
TEMP1 = ONE / TEMP1
*
KILL = LSVEC
LSVEC = RSVEC
RSVEC = KILL
+ IF ( LSVEC ) N1 = N
*
ROWPIV = .TRUE.
END IF
* Assemble the left singular vector matrix U (M x N).
*
IF ( N .LT. M ) THEN
- CALL DLASET( 'A', M-N, N, ZERO, ZERO, U(NR+1,1), LDU )
+ 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(NR+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,
* At this moment, V contains the right singular vectors of A.
* Next, assemble the left singular vector matrix U (M x N).
*
- IF ( N .LT. M ) THEN
- CALL DLASET( 'A', M-N, N, ZERO, ZERO, U(NR+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(NR+1,N+1),LDU )
+ 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 )
END IF
END IF
*
* Arguments
* =========
*
-* JOBA (input) CHARACTER*1
+* JOBA (input) CHARACTER* 1
* Specifies the structure of A.
* = 'L': The input matrix A is lower triangular;
* = 'U': The input matrix A is upper triangular;
* orthogonal up to CTOL*EPS, EPS=DLAMCH('E').
* It is required that CTOL >= ONE, i.e. it is not
* allowed to force the routine to obtain orthogonality
-* below EPSILON.
+* below EPS.
* On exit :
* WORK(1) = SCALE is the scaling factor such that SCALE*SVA(1:N)
* are the computed singular values of A.
* ..
* .. Local Scalars ..
DOUBLE PRECISION AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG,
- + BIGTHETA, CS, CTOL, EPSILON, LARGE, MXAAPQ,
+ + BIGTHETA, CS, CTOL, EPSLN, LARGE, MXAAPQ,
+ MXSINJ, ROOTBIG, ROOTEPS, ROOTSFMIN, ROOTTOL,
- + SCALE, SFMIN, SMALL, SN, T, TEMP1, THETA,
+ + SKL, SFMIN, SMALL, SN, T, TEMP1, THETA,
+ THSIGN, TOL
INTEGER BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1,
+ ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, N2, N34,
* ... and the machine dependent parameters are
*[!] (Make sure that DLAMCH() works properly on the target machine.)
*
- EPSILON = DLAMCH( 'Epsilon' )
- ROOTEPS = DSQRT( EPSILON )
+ EPSLN = DLAMCH( 'Epsilon' )
+ ROOTEPS = DSQRT( EPSLN )
SFMIN = DLAMCH( 'SafeMinimum' )
ROOTSFMIN = DSQRT( SFMIN )
- SMALL = SFMIN / EPSILON
+ SMALL = SFMIN / EPSLN
BIG = DLAMCH( 'Overflow' )
* BIG = ONE / SFMIN
ROOTBIG = ONE / ROOTSFMIN
LARGE = BIG / DSQRT( DBLE( M*N ) )
BIGTHETA = ONE / ROOTEPS
*
- TOL = CTOL*EPSILON
+ TOL = CTOL*EPSLN
ROOTTOL = DSQRT( TOL )
*
- IF( DBLE( M )*EPSILON.GE.ONE ) THEN
+ IF( DBLE( M )*EPSLN.GE.ONE ) THEN
INFO = -5
CALL XERBLA( 'DGESVJ', -INFO )
RETURN
* DSQRT(N)*max_i SVA(i) does not overflow. If INFinite entries
* in A are detected, the procedure returns with INFO=-6.
*
- SCALE = ONE / DSQRT( DBLE( M )*DBLE( N ) )
+ SKL= ONE / DSQRT( DBLE( M )*DBLE( N ) )
NOSCALE = .TRUE.
GOSCALE = .TRUE.
*
* the input matrix is M-by-N lower triangular (trapezoidal)
DO 1874 p = 1, N
AAPP = ZERO
- AAQQ = ZERO
+ AAQQ = ONE
CALL DLASSQ( M-p+1, A( p, p ), 1, AAPP, AAQQ )
IF( AAPP.GT.BIG ) THEN
INFO = -6
SVA( p ) = AAPP*AAQQ
ELSE
NOSCALE = .FALSE.
- SVA( p ) = AAPP*( AAQQ*SCALE )
+ SVA( p ) = AAPP*( AAQQ*SKL)
IF( GOSCALE ) THEN
GOSCALE = .FALSE.
DO 1873 q = 1, p - 1
- SVA( q ) = SVA( q )*SCALE
+ SVA( q ) = SVA( q )*SKL
1873 CONTINUE
END IF
END IF
* the input matrix is M-by-N upper triangular (trapezoidal)
DO 2874 p = 1, N
AAPP = ZERO
- AAQQ = ZERO
+ AAQQ = ONE
CALL DLASSQ( p, A( 1, p ), 1, AAPP, AAQQ )
IF( AAPP.GT.BIG ) THEN
INFO = -6
SVA( p ) = AAPP*AAQQ
ELSE
NOSCALE = .FALSE.
- SVA( p ) = AAPP*( AAQQ*SCALE )
+ SVA( p ) = AAPP*( AAQQ*SKL)
IF( GOSCALE ) THEN
GOSCALE = .FALSE.
DO 2873 q = 1, p - 1
- SVA( q ) = SVA( q )*SCALE
+ SVA( q ) = SVA( q )*SKL
2873 CONTINUE
END IF
END IF
* the input matrix is M-by-N general dense
DO 3874 p = 1, N
AAPP = ZERO
- AAQQ = ZERO
+ AAQQ = ONE
CALL DLASSQ( M, A( 1, p ), 1, AAPP, AAQQ )
IF( AAPP.GT.BIG ) THEN
INFO = -6
SVA( p ) = AAPP*AAQQ
ELSE
NOSCALE = .FALSE.
- SVA( p ) = AAPP*( AAQQ*SCALE )
+ SVA( p ) = AAPP*( AAQQ*SKL)
IF( GOSCALE ) THEN
GOSCALE = .FALSE.
DO 3873 q = 1, p - 1
- SVA( q ) = SVA( q )*SCALE
+ SVA( q ) = SVA( q )*SKL
3873 CONTINUE
END IF
END IF
3874 CONTINUE
END IF
*
- IF( NOSCALE )SCALE = ONE
+ IF( NOSCALE )SKL= ONE
*
* Move the smaller part of the spectrum from the underflow threshold
*(!) Start by determining the position of the nonzero entries of the
* #:) Quick return for one-column matrix
*
IF( N.EQ.1 ) THEN
- IF( LSVEC )CALL DLASCL( 'G', 0, 0, SVA( 1 ), SCALE, M, 1,
+ IF( LSVEC )CALL DLASCL( 'G', 0, 0, SVA( 1 ), SKL, M, 1,
+ A( 1, 1 ), LDA, IERR )
- WORK( 1 ) = ONE / SCALE
+ WORK( 1 ) = ONE / SKL
IF( SVA( 1 ).GE.SFMIN ) THEN
WORK( 2 ) = ONE
ELSE
* Protect small singular values from underflow, and try to
* avoid underflows/overflows in computing Jacobi rotations.
*
- SN = DSQRT( SFMIN / EPSILON )
+ SN = DSQRT( SFMIN / EPSLN )
TEMP1 = DSQRT( BIG / DBLE( N ) )
IF( ( AAPP.LE.SN ) .OR. ( AAQQ.GE.TEMP1 ) .OR.
+ ( ( SN.LE.AAQQ ) .AND. ( AAPP.LE.TEMP1 ) ) ) THEN
IF( TEMP1.NE.ONE ) THEN
CALL DLASCL( 'G', 0, 0, ONE, TEMP1, N, 1, SVA, N, IERR )
END IF
- SCALE = TEMP1*SCALE
- IF( SCALE.NE.ONE ) THEN
- CALL DLASCL( JOBA, 0, 0, ONE, SCALE, M, N, A, LDA, IERR )
- SCALE = ONE / SCALE
+ SKL= TEMP1*SKL
+ IF( SKL.NE.ONE ) THEN
+ CALL DLASCL( JOBA, 0, 0, ONE, SKL, M, N, A, LDA, IERR )
+ SKL= ONE / SKL
END IF
*
* Row-cyclic Jacobi SVD algorithm with column pivoting
*
CALL DGSVJ0( JOBV, M-N34, N-N34, A( N34+1, N34+1 ), LDA,
+ WORK( N34+1 ), SVA( N34+1 ), MVL,
- + V( N34*q+1, N34+1 ), LDV, EPSILON, SFMIN, TOL,
+ + V( N34*q+1, N34+1 ), LDV, EPSLN, SFMIN, TOL,
+ 2, WORK( N+1 ), LWORK-N, IERR )
*
CALL DGSVJ0( JOBV, M-N2, N34-N2, A( N2+1, N2+1 ), LDA,
+ WORK( N2+1 ), SVA( N2+1 ), MVL,
- + V( N2*q+1, N2+1 ), LDV, EPSILON, SFMIN, TOL, 2,
+ + V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 2,
+ WORK( N+1 ), LWORK-N, IERR )
*
CALL DGSVJ1( JOBV, M-N2, N-N2, N4, A( N2+1, N2+1 ), LDA,
+ WORK( N2+1 ), SVA( N2+1 ), MVL,
- + V( N2*q+1, N2+1 ), LDV, EPSILON, SFMIN, TOL, 1,
+ + V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 1,
+ WORK( N+1 ), LWORK-N, IERR )
*
CALL DGSVJ0( JOBV, M-N4, N2-N4, A( N4+1, N4+1 ), LDA,
+ WORK( N4+1 ), SVA( N4+1 ), MVL,
- + V( N4*q+1, N4+1 ), LDV, EPSILON, SFMIN, TOL, 1,
+ + V( N4*q+1, N4+1 ), LDV, EPSLN, SFMIN, TOL, 1,
+ WORK( N+1 ), LWORK-N, IERR )
*
CALL DGSVJ0( JOBV, M, N4, A, LDA, WORK, SVA, MVL, V, LDV,
- + EPSILON, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,
+ + EPSLN, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,
+ IERR )
*
CALL DGSVJ1( JOBV, M, N2, N4, A, LDA, WORK, SVA, MVL, V,
- + LDV, EPSILON, SFMIN, TOL, 1, WORK( N+1 ),
+ + LDV, EPSLN, SFMIN, TOL, 1, WORK( N+1 ),
+ LWORK-N, IERR )
*
*
*
*
CALL DGSVJ0( JOBV, N4, N4, A, LDA, WORK, SVA, MVL, V, LDV,
- + EPSILON, SFMIN, TOL, 2, WORK( N+1 ), LWORK-N,
+ + EPSLN, SFMIN, TOL, 2, WORK( N+1 ), LWORK-N,
+ IERR )
*
CALL DGSVJ0( JOBV, N2, N4, A( 1, N4+1 ), LDA, WORK( N4+1 ),
+ SVA( N4+1 ), MVL, V( N4*q+1, N4+1 ), LDV,
- + EPSILON, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,
+ + EPSLN, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,
+ IERR )
*
CALL DGSVJ1( JOBV, N2, N2, N4, A, LDA, WORK, SVA, MVL, V,
- + LDV, EPSILON, SFMIN, TOL, 1, WORK( N+1 ),
+ + LDV, EPSLN, SFMIN, TOL, 1, WORK( N+1 ),
+ LWORK-N, IERR )
*
CALL DGSVJ0( JOBV, N2+N4, N4, A( 1, N2+1 ), LDA,
+ WORK( N2+1 ), SVA( N2+1 ), MVL,
- + V( N2*q+1, N2+1 ), LDV, EPSILON, SFMIN, TOL, 1,
+ + V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 1,
+ WORK( N+1 ), LWORK-N, IERR )
END IF
SVA( p ) = DNRM2( M, A( 1, p ), 1 )*WORK( p )
ELSE
TEMP1 = ZERO
- AAPP = ZERO
+ AAPP = ONE
CALL DLASSQ( M, A( 1, p ), 1, TEMP1, AAPP )
SVA( p ) = TEMP1*DSQRT( AAPP )*WORK( p )
END IF
+ FASTR )
SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
+ ONE+T*APOAQ*AAPQ ) )
- AAPP = AAPP*DSQRT( ONE-T*AQOAP*
- + AAPQ )
+ AAPP = AAPP*DSQRT( DMAX1( ZERO,
+ + ONE-T*AQOAP*AAPQ ) )
MXSINJ = DMAX1( MXSINJ, DABS( T ) )
*
ELSE
+ WORK( q )
ELSE
T = ZERO
- AAQQ = ZERO
+ AAQQ = ONE
CALL DLASSQ( M, A( 1, q ), 1, T,
+ AAQQ )
SVA( q ) = T*DSQRT( AAQQ )*WORK( q )
+ WORK( p )
ELSE
T = ZERO
- AAPP = ZERO
+ AAPP = ONE
CALL DLASSQ( M, A( 1, p ), 1, T,
+ AAPP )
AAPP = T*DSQRT( AAPP )*WORK( p )
MXSINJ = DMAX1( MXSINJ, DABS( SN ) )
SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
+ ONE+T*APOAQ*AAPQ ) )
- AAPP = AAPP*DSQRT( ONE-T*AQOAP*
- + AAPQ )
+ AAPP = AAPP*DSQRT( DMAX1( ZERO,
+ + ONE-T*AQOAP*AAPQ ) )
*
APOAQ = WORK( p ) / WORK( q )
AQOAP = WORK( q ) / WORK( p )
+ WORK( q )
ELSE
T = ZERO
- AAQQ = ZERO
+ AAQQ = ONE
CALL DLASSQ( M, A( 1, q ), 1, T,
+ AAQQ )
SVA( q ) = T*DSQRT( AAQQ )*WORK( q )
+ WORK( p )
ELSE
T = ZERO
- AAPP = ZERO
+ AAPP = ONE
CALL DLASSQ( M, A( 1, p ), 1, T,
+ AAPP )
AAPP = T*DSQRT( AAPP )*WORK( p )
SVA( N ) = DNRM2( M, A( 1, N ), 1 )*WORK( N )
ELSE
T = ZERO
- AAPP = ZERO
+ AAPP = ONE
CALL DLASSQ( M, A( 1, N ), 1, T, AAPP )
SVA( N ) = T*DSQRT( AAPP )*WORK( N )
END IF
END IF
IF( SVA( p ).NE.ZERO ) THEN
N4 = N4 + 1
- IF( SVA( p )*SCALE.GT.SFMIN )N2 = N2 + 1
+ IF( SVA( p )*SKL.GT.SFMIN )N2 = N2 + 1
END IF
5991 CONTINUE
IF( SVA( N ).NE.ZERO ) THEN
N4 = N4 + 1
- IF( SVA( N )*SCALE.GT.SFMIN )N2 = N2 + 1
+ IF( SVA( N )*SKL.GT.SFMIN )N2 = N2 + 1
END IF
*
* Normalize the left singular vectors.
END IF
*
* Undo scaling, if necessary (and possible).
- IF( ( ( SCALE.GT.ONE ) .AND. ( SVA( 1 ).LT.( BIG /
- + SCALE ) ) ) .OR. ( ( SCALE.LT.ONE ) .AND. ( SVA( N2 ).GT.
- + ( SFMIN / SCALE ) ) ) ) THEN
+ IF( ( ( SKL.GT.ONE ) .AND. ( SVA( 1 ).LT.( BIG /
+ + SKL) ) ) .OR. ( ( SKL.LT.ONE ) .AND. ( SVA( N2 ).GT.
+ + ( SFMIN / SKL) ) ) ) THEN
DO 2400 p = 1, N
- SVA( p ) = SCALE*SVA( p )
+ SVA( p ) = SKL*SVA( p )
2400 CONTINUE
- SCALE = ONE
+ SKL= ONE
END IF
*
- WORK( 1 ) = SCALE
-* The singular values of A are SCALE*SVA(1:N). If SCALE.NE.ONE
+ WORK( 1 ) = SKL
+* The singular values of A are SKL*SVA(1:N). If SKL.NE.ONE
* then some of the singular values may overflow or underflow and
* the spectrum is given in this factored representation.
*
INFO = -3
ELSE IF( LDA.LT.M ) THEN
INFO = -5
- ELSE IF( MV.LT.0 ) THEN
+ ELSE IF( ( RSVEC.OR.APPLV ) .AND. ( MV.LT.0 ) ) THEN
INFO = -8
- ELSE IF( LDV.LT.M ) THEN
+ ELSE IF( ( RSVEC.AND.( LDV.LT.N ) ).OR.
+ & ( APPLV.AND.( LDV.LT.MV ) ) ) THEN
INFO = -10
ELSE IF( TOL.LE.EPS ) THEN
INFO = -13
SVA( p ) = DNRM2( M, A( 1, p ), 1 )*D( p )
ELSE
TEMP1 = ZERO
- AAPP = ZERO
+ AAPP = ONE
CALL DLASSQ( M, A( 1, p ), 1, TEMP1, AAPP )
SVA( p ) = TEMP1*DSQRT( AAPP )*D( p )
END IF
+ FASTR )
SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
+ ONE+T*APOAQ*AAPQ ) )
- AAPP = AAPP*DSQRT( ONE-T*AQOAP*
- + AAPQ )
+ AAPP = AAPP*DSQRT( DMAX1( ZERO,
+ + ONE-T*AQOAP*AAPQ ) )
MXSINJ = DMAX1( MXSINJ, DABS( T ) )
*
ELSE
+ D( q )
ELSE
T = ZERO
- AAQQ = ZERO
+ AAQQ = ONE
CALL DLASSQ( M, A( 1, q ), 1, T,
+ AAQQ )
SVA( q ) = T*DSQRT( AAQQ )*D( q )
+ D( p )
ELSE
T = ZERO
- AAPP = ZERO
+ AAPP = ONE
CALL DLASSQ( M, A( 1, p ), 1, T,
+ AAPP )
AAPP = T*DSQRT( AAPP )*D( p )
MXSINJ = DMAX1( MXSINJ, DABS( SN ) )
SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
+ ONE+T*APOAQ*AAPQ ) )
- AAPP = AAPP*DSQRT( ONE-T*AQOAP*
- + AAPQ )
+ AAPP = AAPP*DSQRT( DMAX1( ZERO,
+ + ONE-T*AQOAP*AAPQ ) )
*
APOAQ = D( p ) / D( q )
AQOAP = D( q ) / D( p )
+ D( q )
ELSE
T = ZERO
- AAQQ = ZERO
+ AAQQ = ONE
CALL DLASSQ( M, A( 1, q ), 1, T,
+ AAQQ )
SVA( q ) = T*DSQRT( AAQQ )*D( q )
+ D( p )
ELSE
T = ZERO
- AAPP = ZERO
+ AAPP = ONE
CALL DLASSQ( M, A( 1, p ), 1, T,
+ AAPP )
AAPP = T*DSQRT( AAPP )*D( p )
SVA( N ) = DNRM2( M, A( 1, N ), 1 )*D( N )
ELSE
T = ZERO
- AAPP = ZERO
+ AAPP = ONE
CALL DLASSQ( M, A( 1, N ), 1, T, AAPP )
SVA( N ) = T*DSQRT( AAPP )*D( N )
END IF
INFO = -4
ELSE IF( LDA.LT.M ) THEN
INFO = -6
- ELSE IF( MV.LT.0 ) THEN
+ ELSE IF( ( RSVEC.OR.APPLV ) .AND. ( MV.LT.0 ) ) THEN
INFO = -9
- ELSE IF( LDV.LT.M ) THEN
+ ELSE IF( ( RSVEC.AND.( LDV.LT.N ) ).OR.
+ & ( APPLV.AND.( LDV.LT.MV ) ) ) THEN
INFO = -11
ELSE IF( TOL.LE.EPS ) THEN
INFO = -14
MXSINJ = DMAX1( MXSINJ, DABS( SN ) )
SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
+ ONE+T*APOAQ*AAPQ ) )
- AAPP = AAPP*DSQRT( ONE-T*AQOAP*
- + AAPQ )
+ AAPP = AAPP*DSQRT( DMAX1( ZERO,
+ + ONE-T*AQOAP*AAPQ ) )
APOAQ = D( p ) / D( q )
AQOAP = D( q ) / D( p )
+ D( q )
ELSE
T = ZERO
- AAQQ = ZERO
+ AAQQ = ONE
CALL DLASSQ( M, A( 1, q ), 1, T,
+ AAQQ )
SVA( q ) = T*DSQRT( AAQQ )*D( q )
+ D( p )
ELSE
T = ZERO
- AAPP = ZERO
+ AAPP = ONE
CALL DLASSQ( M, A( 1, p ), 1, T,
+ AAPP )
AAPP = T*DSQRT( AAPP )*D( p )
SVA( N ) = DNRM2( M, A( 1, N ), 1 )*D( N )
ELSE
T = ZERO
- AAPP = ZERO
+ AAPP = ONE
CALL DLASSQ( M, A( 1, N ), 1, T, AAPP )
SVA( N ) = T*DSQRT( AAPP )*D( N )
END IF
GOSCAL = .TRUE.
DO 1874 p = 1, N
AAPP = ZERO
- AAQQ = ZERO
+ AAQQ = ONE
CALL SLASSQ( M, A(1,p), 1, AAPP, AAQQ )
IF ( AAPP .GT. BIG ) THEN
INFO = - 9
IF ( L2TRAN ) THEN
DO 1950 p = 1, M
XSC = ZERO
- TEMP1 = ZERO
+ TEMP1 = ONE
CALL SLASSQ( N, A(p,1), LDA, XSC, TEMP1 )
* SLASSQ gets both the ell_2 and the ell_infinity norm
* in one pass through the vector
IF ( L2TRAN ) THEN
*
XSC = ZERO
- TEMP1 = ZERO
+ TEMP1 = ONE
CALL SLASSQ( N, SVA, 1, XSC, TEMP1 )
TEMP1 = ONE / TEMP1
*
KILL = LSVEC
LSVEC = RSVEC
RSVEC = KILL
+ IF ( LSVEC ) N1 = N
*
ROWPIV = .TRUE.
END IF
* Assemble the left singular vector matrix U (M x N).
*
IF ( N .LT. M ) THEN
- CALL SLASET( 'A', M-N, N, ZERO, ZERO, U(NR+1,1), LDU )
+ 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(NR+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,
* At this moment, V contains the right singular vectors of A.
* Next, assemble the left singular vector matrix U (M x N).
*
- IF ( N .LT. M ) THEN
- CALL SLASET( 'A', M-N, N, ZERO, ZERO, U(NR+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(NR+1,N+1),LDU )
+ 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 )
END IF
END IF
*
* Arguments
* =========
*
-* JOBA (input) CHARACTER*1
+* JOBA (input) CHARACTER* 1
* Specifies the structure of A.
* = 'L': The input matrix A is lower triangular;
* = 'U': The input matrix A is upper triangular;
* ..
* .. Local Scalars ..
REAL AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG,
- + BIGTHETA, CS, CTOL, EPSILON, LARGE, MXAAPQ,
+ + BIGTHETA, CS, CTOL, EPSLN, LARGE, MXAAPQ,
+ MXSINJ, ROOTBIG, ROOTEPS, ROOTSFMIN, ROOTTOL,
- + SCALE, SFMIN, SMALL, SN, T, TEMP1, THETA,
+ + SKL, SFMIN, SMALL, SN, T, TEMP1, THETA,
+ THSIGN, TOL
INTEGER BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1,
+ ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, N2, N34,
* ... and the machine dependent parameters are
*[!] (Make sure that SLAMCH() works properly on the target machine.)
*
- EPSILON = SLAMCH( 'Epsilon' )
- ROOTEPS = SQRT( EPSILON )
+ EPSLN = SLAMCH( 'Epsilon' )
+ ROOTEPS = SQRT( EPSLN )
SFMIN = SLAMCH( 'SafeMinimum' )
ROOTSFMIN = SQRT( SFMIN )
- SMALL = SFMIN / EPSILON
+ SMALL = SFMIN / EPSLN
BIG = SLAMCH( 'Overflow' )
ROOTBIG = ONE / ROOTSFMIN
LARGE = BIG / SQRT( FLOAT( M*N ) )
BIGTHETA = ONE / ROOTEPS
*
- TOL = CTOL*EPSILON
+ TOL = CTOL*EPSLN
ROOTTOL = SQRT( TOL )
*
- IF( FLOAT( M )*EPSILON.GE.ONE ) THEN
+ IF( FLOAT( M )*EPSLN.GE.ONE ) THEN
INFO = -5
CALL XERBLA( 'SGESVJ', -INFO )
RETURN
* SQRT(N)*max_i SVA(i) does not overflow. If INFinite entries
* in A are detected, the procedure returns with INFO=-6.
*
- SCALE = ONE / SQRT( FLOAT( M )*FLOAT( N ) )
+ SKL = ONE / SQRT( FLOAT( M )*FLOAT( N ) )
NOSCALE = .TRUE.
GOSCALE = .TRUE.
*
* the input matrix is M-by-N lower triangular (trapezoidal)
DO 1874 p = 1, N
AAPP = ZERO
- AAQQ = ZERO
+ AAQQ = ONE
CALL SLASSQ( M-p+1, A( p, p ), 1, AAPP, AAQQ )
IF( AAPP.GT.BIG ) THEN
INFO = -6
SVA( p ) = AAPP*AAQQ
ELSE
NOSCALE = .FALSE.
- SVA( p ) = AAPP*( AAQQ*SCALE )
+ SVA( p ) = AAPP*( AAQQ*SKL )
IF( GOSCALE ) THEN
GOSCALE = .FALSE.
DO 1873 q = 1, p - 1
- SVA( q ) = SVA( q )*SCALE
+ SVA( q ) = SVA( q )*SKL
1873 CONTINUE
END IF
END IF
* the input matrix is M-by-N upper triangular (trapezoidal)
DO 2874 p = 1, N
AAPP = ZERO
- AAQQ = ZERO
+ AAQQ = ONE
CALL SLASSQ( p, A( 1, p ), 1, AAPP, AAQQ )
IF( AAPP.GT.BIG ) THEN
INFO = -6
SVA( p ) = AAPP*AAQQ
ELSE
NOSCALE = .FALSE.
- SVA( p ) = AAPP*( AAQQ*SCALE )
+ SVA( p ) = AAPP*( AAQQ*SKL )
IF( GOSCALE ) THEN
GOSCALE = .FALSE.
DO 2873 q = 1, p - 1
- SVA( q ) = SVA( q )*SCALE
+ SVA( q ) = SVA( q )*SKL
2873 CONTINUE
END IF
END IF
* the input matrix is M-by-N general dense
DO 3874 p = 1, N
AAPP = ZERO
- AAQQ = ZERO
+ AAQQ = ONE
CALL SLASSQ( M, A( 1, p ), 1, AAPP, AAQQ )
IF( AAPP.GT.BIG ) THEN
INFO = -6
SVA( p ) = AAPP*AAQQ
ELSE
NOSCALE = .FALSE.
- SVA( p ) = AAPP*( AAQQ*SCALE )
+ SVA( p ) = AAPP*( AAQQ*SKL )
IF( GOSCALE ) THEN
GOSCALE = .FALSE.
DO 3873 q = 1, p - 1
- SVA( q ) = SVA( q )*SCALE
+ SVA( q ) = SVA( q )*SKL
3873 CONTINUE
END IF
END IF
3874 CONTINUE
END IF
*
- IF( NOSCALE )SCALE = ONE
+ IF( NOSCALE )SKL = ONE
*
* Move the smaller part of the spectrum from the underflow threshold
*(!) Start by determining the position of the nonzero entries of the
* #:) Quick return for one-column matrix
*
IF( N.EQ.1 ) THEN
- IF( LSVEC )CALL SLASCL( 'G', 0, 0, SVA( 1 ), SCALE, M, 1,
+ IF( LSVEC )CALL SLASCL( 'G', 0, 0, SVA( 1 ), SKL, M, 1,
+ A( 1, 1 ), LDA, IERR )
- WORK( 1 ) = ONE / SCALE
+ WORK( 1 ) = ONE / SKL
IF( SVA( 1 ).GE.SFMIN ) THEN
WORK( 2 ) = ONE
ELSE
* Protect small singular values from underflow, and try to
* avoid underflows/overflows in computing Jacobi rotations.
*
- SN = SQRT( SFMIN / EPSILON )
+ SN = SQRT( SFMIN / EPSLN )
TEMP1 = SQRT( BIG / FLOAT( N ) )
IF( ( AAPP.LE.SN ) .OR. ( AAQQ.GE.TEMP1 ) .OR.
+ ( ( SN.LE.AAQQ ) .AND. ( AAPP.LE.TEMP1 ) ) ) THEN
IF( TEMP1.NE.ONE ) THEN
CALL SLASCL( 'G', 0, 0, ONE, TEMP1, N, 1, SVA, N, IERR )
END IF
- SCALE = TEMP1*SCALE
- IF( SCALE.NE.ONE ) THEN
- CALL SLASCL( JOBA, 0, 0, ONE, SCALE, M, N, A, LDA, IERR )
- SCALE = ONE / SCALE
+ SKL = TEMP1*SKL
+ IF( SKL.NE.ONE ) THEN
+ CALL SLASCL( JOBA, 0, 0, ONE, SKL, M, N, A, LDA, IERR )
+ SKL = ONE / SKL
END IF
*
* Row-cyclic Jacobi SVD algorithm with column pivoting
*
CALL SGSVJ0( JOBV, M-N34, N-N34, A( N34+1, N34+1 ), LDA,
+ WORK( N34+1 ), SVA( N34+1 ), MVL,
- + V( N34*q+1, N34+1 ), LDV, EPSILON, SFMIN, TOL,
+ + V( N34*q+1, N34+1 ), LDV, EPSLN, SFMIN, TOL,
+ 2, WORK( N+1 ), LWORK-N, IERR )
*
CALL SGSVJ0( JOBV, M-N2, N34-N2, A( N2+1, N2+1 ), LDA,
+ WORK( N2+1 ), SVA( N2+1 ), MVL,
- + V( N2*q+1, N2+1 ), LDV, EPSILON, SFMIN, TOL, 2,
+ + V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 2,
+ WORK( N+1 ), LWORK-N, IERR )
*
CALL SGSVJ1( JOBV, M-N2, N-N2, N4, A( N2+1, N2+1 ), LDA,
+ WORK( N2+1 ), SVA( N2+1 ), MVL,
- + V( N2*q+1, N2+1 ), LDV, EPSILON, SFMIN, TOL, 1,
+ + V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 1,
+ WORK( N+1 ), LWORK-N, IERR )
*
CALL SGSVJ0( JOBV, M-N4, N2-N4, A( N4+1, N4+1 ), LDA,
+ WORK( N4+1 ), SVA( N4+1 ), MVL,
- + V( N4*q+1, N4+1 ), LDV, EPSILON, SFMIN, TOL, 1,
+ + V( N4*q+1, N4+1 ), LDV, EPSLN, SFMIN, TOL, 1,
+ WORK( N+1 ), LWORK-N, IERR )
*
CALL SGSVJ0( JOBV, M, N4, A, LDA, WORK, SVA, MVL, V, LDV,
- + EPSILON, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,
+ + EPSLN, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,
+ IERR )
*
CALL SGSVJ1( JOBV, M, N2, N4, A, LDA, WORK, SVA, MVL, V,
- + LDV, EPSILON, SFMIN, TOL, 1, WORK( N+1 ),
+ + LDV, EPSLN, SFMIN, TOL, 1, WORK( N+1 ),
+ LWORK-N, IERR )
*
*
*
*
CALL SGSVJ0( JOBV, N4, N4, A, LDA, WORK, SVA, MVL, V, LDV,
- + EPSILON, SFMIN, TOL, 2, WORK( N+1 ), LWORK-N,
+ + EPSLN, SFMIN, TOL, 2, WORK( N+1 ), LWORK-N,
+ IERR )
*
CALL SGSVJ0( JOBV, N2, N4, A( 1, N4+1 ), LDA, WORK( N4+1 ),
+ SVA( N4+1 ), MVL, V( N4*q+1, N4+1 ), LDV,
- + EPSILON, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,
+ + EPSLN, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,
+ IERR )
*
CALL SGSVJ1( JOBV, N2, N2, N4, A, LDA, WORK, SVA, MVL, V,
- + LDV, EPSILON, SFMIN, TOL, 1, WORK( N+1 ),
+ + LDV, EPSLN, SFMIN, TOL, 1, WORK( N+1 ),
+ LWORK-N, IERR )
*
CALL SGSVJ0( JOBV, N2+N4, N4, A( 1, N2+1 ), LDA,
+ WORK( N2+1 ), SVA( N2+1 ), MVL,
- + V( N2*q+1, N2+1 ), LDV, EPSILON, SFMIN, TOL, 1,
+ + V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 1,
+ WORK( N+1 ), LWORK-N, IERR )
END IF
SVA( p ) = SNRM2( M, A( 1, p ), 1 )*WORK( p )
ELSE
TEMP1 = ZERO
- AAPP = ZERO
+ AAPP = ONE
CALL SLASSQ( M, A( 1, p ), 1, TEMP1, AAPP )
SVA( p ) = TEMP1*SQRT( AAPP )*WORK( p )
END IF
+ FASTR )
SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
+ ONE+T*APOAQ*AAPQ ) )
- AAPP = AAPP*SQRT( ONE-T*AQOAP*AAPQ )
+ AAPP = AAPP*SQRT( AMAX1( ZERO,
+ + ONE-T*AQOAP*AAPQ ) )
MXSINJ = AMAX1( MXSINJ, ABS( T ) )
*
ELSE
+ WORK( q )
ELSE
T = ZERO
- AAQQ = ZERO
+ AAQQ = ONE
CALL SLASSQ( M, A( 1, q ), 1, T,
+ AAQQ )
SVA( q ) = T*SQRT( AAQQ )*WORK( q )
+ WORK( p )
ELSE
T = ZERO
- AAPP = ZERO
+ AAPP = ONE
CALL SLASSQ( M, A( 1, p ), 1, T,
+ AAPP )
AAPP = T*SQRT( AAPP )*WORK( p )
MXSINJ = AMAX1( MXSINJ, ABS( SN ) )
SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
+ ONE+T*APOAQ*AAPQ ) )
- AAPP = AAPP*SQRT( ONE-T*AQOAP*AAPQ )
+ AAPP = AAPP*SQRT( AMAX1( ZERO,
+ + ONE-T*AQOAP*AAPQ ) )
*
APOAQ = WORK( p ) / WORK( q )
AQOAP = WORK( q ) / WORK( p )
+ WORK( q )
ELSE
T = ZERO
- AAQQ = ZERO
+ AAQQ = ONE
CALL SLASSQ( M, A( 1, q ), 1, T,
+ AAQQ )
SVA( q ) = T*SQRT( AAQQ )*WORK( q )
+ WORK( p )
ELSE
T = ZERO
- AAPP = ZERO
+ AAPP = ONE
CALL SLASSQ( M, A( 1, p ), 1, T,
+ AAPP )
AAPP = T*SQRT( AAPP )*WORK( p )
SVA( N ) = SNRM2( M, A( 1, N ), 1 )*WORK( N )
ELSE
T = ZERO
- AAPP = ZERO
+ AAPP = ONE
CALL SLASSQ( M, A( 1, N ), 1, T, AAPP )
SVA( N ) = T*SQRT( AAPP )*WORK( N )
END IF
END IF
IF( SVA( p ).NE.ZERO ) THEN
N4 = N4 + 1
- IF( SVA( p )*SCALE.GT.SFMIN )N2 = N2 + 1
+ IF( SVA( p )*SKL.GT.SFMIN )N2 = N2 + 1
END IF
5991 CONTINUE
IF( SVA( N ).NE.ZERO ) THEN
N4 = N4 + 1
- IF( SVA( N )*SCALE.GT.SFMIN )N2 = N2 + 1
+ IF( SVA( N )*SKL.GT.SFMIN )N2 = N2 + 1
END IF
*
* Normalize the left singular vectors.
END IF
*
* Undo scaling, if necessary (and possible).
- IF( ( ( SCALE.GT.ONE ) .AND. ( SVA( 1 ).LT.( BIG /
- + SCALE ) ) ) .OR. ( ( SCALE.LT.ONE ) .AND. ( SVA( N2 ).GT.
- + ( SFMIN / SCALE ) ) ) ) THEN
+ IF( ( ( SKL.GT.ONE ) .AND. ( SVA( 1 ).LT.( BIG /
+ + SKL ) ) ) .OR. ( ( SKL.LT.ONE ) .AND. ( SVA( N2 ).GT.
+ + ( SFMIN / SKL ) ) ) ) THEN
DO 2400 p = 1, N
- SVA( p ) = SCALE*SVA( p )
+ SVA( p ) = SKL*SVA( p )
2400 CONTINUE
- SCALE = ONE
+ SKL = ONE
END IF
*
- WORK( 1 ) = SCALE
-* The singular values of A are SCALE*SVA(1:N). If SCALE.NE.ONE
+ WORK( 1 ) = SKL
+* The singular values of A are SKL*SVA(1:N). If SKL.NE.ONE
* then some of the singular values may overflow or underflow and
* the spectrum is given in this factored representation.
*
INFO = -3
ELSE IF( LDA.LT.M ) THEN
INFO = -5
- ELSE IF( MV.LT.0 ) THEN
+ ELSE IF( ( RSVEC.OR.APPLV ) .AND. ( MV.LT.0 ) ) THEN
INFO = -8
- ELSE IF( LDV.LT.M ) THEN
+ ELSE IF( ( RSVEC.AND.( LDV.LT.N ) ).OR.
+ & ( APPLV.AND.( LDV.LT.MV ) ) ) THEN
INFO = -10
ELSE IF( TOL.LE.EPS ) THEN
INFO = -13
SVA( p ) = SNRM2( M, A( 1, p ), 1 )*D( p )
ELSE
TEMP1 = ZERO
- AAPP = ZERO
+ AAPP = ONE
CALL SLASSQ( M, A( 1, p ), 1, TEMP1, AAPP )
SVA( p ) = TEMP1*SQRT( AAPP )*D( p )
END IF
+ FASTR )
SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
+ ONE+T*APOAQ*AAPQ ) )
- AAPP = AAPP*SQRT( ONE-T*AQOAP*AAPQ )
+ AAPP = AAPP*SQRT( AMAX1( ZERO,
+ + ONE-T*AQOAP*AAPQ ) )
MXSINJ = AMAX1( MXSINJ, ABS( T ) )
*
ELSE
+ D( q )
ELSE
T = ZERO
- AAQQ = ZERO
+ AAQQ = ONE
CALL SLASSQ( M, A( 1, q ), 1, T,
+ AAQQ )
SVA( q ) = T*SQRT( AAQQ )*D( q )
+ D( p )
ELSE
T = ZERO
- AAPP = ZERO
+ AAPP = ONE
CALL SLASSQ( M, A( 1, p ), 1, T,
+ AAPP )
AAPP = T*SQRT( AAPP )*D( p )
MXSINJ = AMAX1( MXSINJ, ABS( SN ) )
SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
+ ONE+T*APOAQ*AAPQ ) )
- AAPP = AAPP*SQRT( ONE-T*AQOAP*AAPQ )
+ AAPP = AAPP*SQRT( AMAX1( ZERO,
+ + ONE-T*AQOAP*AAPQ ) )
*
APOAQ = D( p ) / D( q )
AQOAP = D( q ) / D( p )
+ D( q )
ELSE
T = ZERO
- AAQQ = ZERO
+ AAQQ = ONE
CALL SLASSQ( M, A( 1, q ), 1, T,
+ AAQQ )
SVA( q ) = T*SQRT( AAQQ )*D( q )
+ D( p )
ELSE
T = ZERO
- AAPP = ZERO
+ AAPP = ONE
CALL SLASSQ( M, A( 1, p ), 1, T,
+ AAPP )
AAPP = T*SQRT( AAPP )*D( p )
SVA( N ) = SNRM2( M, A( 1, N ), 1 )*D( N )
ELSE
T = ZERO
- AAPP = ZERO
+ AAPP = ONE
CALL SLASSQ( M, A( 1, N ), 1, T, AAPP )
SVA( N ) = T*SQRT( AAPP )*D( N )
END IF
INFO = -4
ELSE IF( LDA.LT.M ) THEN
INFO = -6
- ELSE IF( MV.LT.0 ) THEN
+ ELSE IF( ( RSVEC.OR.APPLV ) .AND. ( MV.LT.0 ) ) THEN
INFO = -9
- ELSE IF( LDV.LT.M ) THEN
+ ELSE IF( ( RSVEC.AND.( LDV.LT.N ) ).OR.
+ & ( APPLV.AND.( LDV.LT.MV ) ) ) THEN
INFO = -11
ELSE IF( TOL.LE.EPS ) THEN
INFO = -14
MXSINJ = AMAX1( MXSINJ, ABS( SN ) )
SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
+ ONE+T*APOAQ*AAPQ ) )
- AAPP = AAPP*SQRT( ONE-T*AQOAP*AAPQ )
+ AAPP = AAPP*SQRT( AMAX1( ZERO,
+ + ONE-T*AQOAP*AAPQ ) )
APOAQ = D( p ) / D( q )
AQOAP = D( q ) / D( p )
+ D( q )
ELSE
T = ZERO
- AAQQ = ZERO
+ AAQQ = ONE
CALL SLASSQ( M, A( 1, q ), 1, T,
+ AAQQ )
SVA( q ) = T*SQRT( AAQQ )*D( q )
+ D( p )
ELSE
T = ZERO
- AAPP = ZERO
+ AAPP = ONE
CALL SLASSQ( M, A( 1, p ), 1, T,
+ AAPP )
AAPP = T*SQRT( AAPP )*D( p )
SVA( N ) = SNRM2( M, A( 1, N ), 1 )*D( N )
ELSE
T = ZERO
- AAPP = ZERO
+ AAPP = ONE
CALL SLASSQ( M, A( 1, N ), 1, T, AAPP )
SVA( N ) = T*SQRT( AAPP )*D( N )
END IF