version 1.1, 2010/08/07 13:21:03
|
version 1.5, 2010/12/21 13:53:26
|
Line 1
|
Line 1
|
SUBROUTINE DGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, |
SUBROUTINE DGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, |
+ LDV, WORK, LWORK, INFO ) |
+ LDV, WORK, LWORK, INFO ) |
* |
* |
* -- LAPACK routine (version 3.2.2) -- |
* -- LAPACK routine (version 3.3.0) -- |
* |
* |
* -- Contributed by Zlatko Drmac of the University of Zagreb and -- |
* -- Contributed by Zlatko Drmac of the University of Zagreb and -- |
* -- Kresimir Veselic of the Fernuniversitaet Hagen -- |
* -- Kresimir Veselic of the Fernuniversitaet Hagen -- |
* -- June 2010 -- |
* November 2010 |
* |
* |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
Line 222
|
Line 222
|
* orthogonal up to CTOL*EPS, EPS=DLAMCH('E'). |
* orthogonal up to CTOL*EPS, EPS=DLAMCH('E'). |
* It is required that CTOL >= ONE, i.e. it is not |
* It is required that CTOL >= ONE, i.e. it is not |
* allowed to force the routine to obtain orthogonality |
* allowed to force the routine to obtain orthogonality |
* below EPSILON. |
* below EPS. |
* On exit : |
* On exit : |
* WORK(1) = SCALE is the scaling factor such that SCALE*SVA(1:N) |
* WORK(1) = SCALE is the scaling factor such that SCALE*SVA(1:N) |
* are the computed singular values of A. |
* are the computed singular values of A. |
Line 262
|
Line 262
|
* .. |
* .. |
* .. Local Scalars .. |
* .. Local Scalars .. |
DOUBLE PRECISION AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG, |
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, |
+ MXSINJ, ROOTBIG, ROOTEPS, ROOTSFMIN, ROOTTOL, |
+ SCALE, SFMIN, SMALL, SN, T, TEMP1, THETA, |
+ SKL, SFMIN, SMALL, SN, T, TEMP1, THETA, |
+ THSIGN, TOL |
+ THSIGN, TOL |
INTEGER BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1, |
INTEGER BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1, |
+ ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, N2, N34, |
+ ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, N2, N34, |
Line 368
|
Line 368
|
* ... and the machine dependent parameters are |
* ... and the machine dependent parameters are |
*[!] (Make sure that DLAMCH() works properly on the target machine.) |
*[!] (Make sure that DLAMCH() works properly on the target machine.) |
* |
* |
EPSILON = DLAMCH( 'Epsilon' ) |
EPSLN = DLAMCH( 'Epsilon' ) |
ROOTEPS = DSQRT( EPSILON ) |
ROOTEPS = DSQRT( EPSLN ) |
SFMIN = DLAMCH( 'SafeMinimum' ) |
SFMIN = DLAMCH( 'SafeMinimum' ) |
ROOTSFMIN = DSQRT( SFMIN ) |
ROOTSFMIN = DSQRT( SFMIN ) |
SMALL = SFMIN / EPSILON |
SMALL = SFMIN / EPSLN |
BIG = DLAMCH( 'Overflow' ) |
BIG = DLAMCH( 'Overflow' ) |
* BIG = ONE / SFMIN |
* BIG = ONE / SFMIN |
ROOTBIG = ONE / ROOTSFMIN |
ROOTBIG = ONE / ROOTSFMIN |
LARGE = BIG / DSQRT( DBLE( M*N ) ) |
LARGE = BIG / DSQRT( DBLE( M*N ) ) |
BIGTHETA = ONE / ROOTEPS |
BIGTHETA = ONE / ROOTEPS |
* |
* |
TOL = CTOL*EPSILON |
TOL = CTOL*EPSLN |
ROOTTOL = DSQRT( TOL ) |
ROOTTOL = DSQRT( TOL ) |
* |
* |
IF( DBLE( M )*EPSILON.GE.ONE ) THEN |
IF( DBLE( M )*EPSLN.GE.ONE ) THEN |
INFO = -5 |
INFO = -5 |
CALL XERBLA( 'DGESVJ', -INFO ) |
CALL XERBLA( 'DGESVJ', -INFO ) |
RETURN |
RETURN |
Line 407
|
Line 407
|
* DSQRT(N)*max_i SVA(i) does not overflow. If INFinite entries |
* DSQRT(N)*max_i SVA(i) does not overflow. If INFinite entries |
* in A are detected, the procedure returns with INFO=-6. |
* 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. |
NOSCALE = .TRUE. |
GOSCALE = .TRUE. |
GOSCALE = .TRUE. |
* |
* |
Line 415
|
Line 415
|
* the input matrix is M-by-N lower triangular (trapezoidal) |
* the input matrix is M-by-N lower triangular (trapezoidal) |
DO 1874 p = 1, N |
DO 1874 p = 1, N |
AAPP = ZERO |
AAPP = ZERO |
AAQQ = ZERO |
AAQQ = ONE |
CALL DLASSQ( M-p+1, A( p, p ), 1, AAPP, AAQQ ) |
CALL DLASSQ( M-p+1, A( p, p ), 1, AAPP, AAQQ ) |
IF( AAPP.GT.BIG ) THEN |
IF( AAPP.GT.BIG ) THEN |
INFO = -6 |
INFO = -6 |
Line 427
|
Line 427
|
SVA( p ) = AAPP*AAQQ |
SVA( p ) = AAPP*AAQQ |
ELSE |
ELSE |
NOSCALE = .FALSE. |
NOSCALE = .FALSE. |
SVA( p ) = AAPP*( AAQQ*SCALE ) |
SVA( p ) = AAPP*( AAQQ*SKL) |
IF( GOSCALE ) THEN |
IF( GOSCALE ) THEN |
GOSCALE = .FALSE. |
GOSCALE = .FALSE. |
DO 1873 q = 1, p - 1 |
DO 1873 q = 1, p - 1 |
SVA( q ) = SVA( q )*SCALE |
SVA( q ) = SVA( q )*SKL |
1873 CONTINUE |
1873 CONTINUE |
END IF |
END IF |
END IF |
END IF |
Line 440
|
Line 440
|
* the input matrix is M-by-N upper triangular (trapezoidal) |
* the input matrix is M-by-N upper triangular (trapezoidal) |
DO 2874 p = 1, N |
DO 2874 p = 1, N |
AAPP = ZERO |
AAPP = ZERO |
AAQQ = ZERO |
AAQQ = ONE |
CALL DLASSQ( p, A( 1, p ), 1, AAPP, AAQQ ) |
CALL DLASSQ( p, A( 1, p ), 1, AAPP, AAQQ ) |
IF( AAPP.GT.BIG ) THEN |
IF( AAPP.GT.BIG ) THEN |
INFO = -6 |
INFO = -6 |
Line 452
|
Line 452
|
SVA( p ) = AAPP*AAQQ |
SVA( p ) = AAPP*AAQQ |
ELSE |
ELSE |
NOSCALE = .FALSE. |
NOSCALE = .FALSE. |
SVA( p ) = AAPP*( AAQQ*SCALE ) |
SVA( p ) = AAPP*( AAQQ*SKL) |
IF( GOSCALE ) THEN |
IF( GOSCALE ) THEN |
GOSCALE = .FALSE. |
GOSCALE = .FALSE. |
DO 2873 q = 1, p - 1 |
DO 2873 q = 1, p - 1 |
SVA( q ) = SVA( q )*SCALE |
SVA( q ) = SVA( q )*SKL |
2873 CONTINUE |
2873 CONTINUE |
END IF |
END IF |
END IF |
END IF |
Line 465
|
Line 465
|
* the input matrix is M-by-N general dense |
* the input matrix is M-by-N general dense |
DO 3874 p = 1, N |
DO 3874 p = 1, N |
AAPP = ZERO |
AAPP = ZERO |
AAQQ = ZERO |
AAQQ = ONE |
CALL DLASSQ( M, A( 1, p ), 1, AAPP, AAQQ ) |
CALL DLASSQ( M, A( 1, p ), 1, AAPP, AAQQ ) |
IF( AAPP.GT.BIG ) THEN |
IF( AAPP.GT.BIG ) THEN |
INFO = -6 |
INFO = -6 |
Line 477
|
Line 477
|
SVA( p ) = AAPP*AAQQ |
SVA( p ) = AAPP*AAQQ |
ELSE |
ELSE |
NOSCALE = .FALSE. |
NOSCALE = .FALSE. |
SVA( p ) = AAPP*( AAQQ*SCALE ) |
SVA( p ) = AAPP*( AAQQ*SKL) |
IF( GOSCALE ) THEN |
IF( GOSCALE ) THEN |
GOSCALE = .FALSE. |
GOSCALE = .FALSE. |
DO 3873 q = 1, p - 1 |
DO 3873 q = 1, p - 1 |
SVA( q ) = SVA( q )*SCALE |
SVA( q ) = SVA( q )*SKL |
3873 CONTINUE |
3873 CONTINUE |
END IF |
END IF |
END IF |
END IF |
3874 CONTINUE |
3874 CONTINUE |
END IF |
END IF |
* |
* |
IF( NOSCALE )SCALE = ONE |
IF( NOSCALE )SKL= ONE |
* |
* |
* Move the smaller part of the spectrum from the underflow threshold |
* Move the smaller part of the spectrum from the underflow threshold |
*(!) Start by determining the position of the nonzero entries of the |
*(!) Start by determining the position of the nonzero entries of the |
Line 517
|
Line 517
|
* #:) Quick return for one-column matrix |
* #:) Quick return for one-column matrix |
* |
* |
IF( N.EQ.1 ) THEN |
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 ) |
+ A( 1, 1 ), LDA, IERR ) |
WORK( 1 ) = ONE / SCALE |
WORK( 1 ) = ONE / SKL |
IF( SVA( 1 ).GE.SFMIN ) THEN |
IF( SVA( 1 ).GE.SFMIN ) THEN |
WORK( 2 ) = ONE |
WORK( 2 ) = ONE |
ELSE |
ELSE |
Line 535
|
Line 535
|
* Protect small singular values from underflow, and try to |
* Protect small singular values from underflow, and try to |
* avoid underflows/overflows in computing Jacobi rotations. |
* avoid underflows/overflows in computing Jacobi rotations. |
* |
* |
SN = DSQRT( SFMIN / EPSILON ) |
SN = DSQRT( SFMIN / EPSLN ) |
TEMP1 = DSQRT( BIG / DBLE( N ) ) |
TEMP1 = DSQRT( BIG / DBLE( N ) ) |
IF( ( AAPP.LE.SN ) .OR. ( AAQQ.GE.TEMP1 ) .OR. |
IF( ( AAPP.LE.SN ) .OR. ( AAQQ.GE.TEMP1 ) .OR. |
+ ( ( SN.LE.AAQQ ) .AND. ( AAPP.LE.TEMP1 ) ) ) THEN |
+ ( ( SN.LE.AAQQ ) .AND. ( AAPP.LE.TEMP1 ) ) ) THEN |
Line 563
|
Line 563
|
IF( TEMP1.NE.ONE ) THEN |
IF( TEMP1.NE.ONE ) THEN |
CALL DLASCL( 'G', 0, 0, ONE, TEMP1, N, 1, SVA, N, IERR ) |
CALL DLASCL( 'G', 0, 0, ONE, TEMP1, N, 1, SVA, N, IERR ) |
END IF |
END IF |
SCALE = TEMP1*SCALE |
SKL= TEMP1*SKL |
IF( SCALE.NE.ONE ) THEN |
IF( SKL.NE.ONE ) THEN |
CALL DLASCL( JOBA, 0, 0, ONE, SCALE, M, N, A, LDA, IERR ) |
CALL DLASCL( JOBA, 0, 0, ONE, SKL, M, N, A, LDA, IERR ) |
SCALE = ONE / SCALE |
SKL= ONE / SKL |
END IF |
END IF |
* |
* |
* Row-cyclic Jacobi SVD algorithm with column pivoting |
* Row-cyclic Jacobi SVD algorithm with column pivoting |
Line 639
|
Line 639
|
* |
* |
CALL DGSVJ0( JOBV, M-N34, N-N34, A( N34+1, N34+1 ), LDA, |
CALL DGSVJ0( JOBV, M-N34, N-N34, A( N34+1, N34+1 ), LDA, |
+ WORK( N34+1 ), SVA( N34+1 ), MVL, |
+ 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 ) |
+ 2, WORK( N+1 ), LWORK-N, IERR ) |
* |
* |
CALL DGSVJ0( JOBV, M-N2, N34-N2, A( N2+1, N2+1 ), LDA, |
CALL DGSVJ0( JOBV, M-N2, N34-N2, A( N2+1, N2+1 ), LDA, |
+ WORK( N2+1 ), SVA( N2+1 ), MVL, |
+ 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 ) |
+ WORK( N+1 ), LWORK-N, IERR ) |
* |
* |
CALL DGSVJ1( JOBV, M-N2, N-N2, N4, A( N2+1, N2+1 ), LDA, |
CALL DGSVJ1( JOBV, M-N2, N-N2, N4, A( N2+1, N2+1 ), LDA, |
+ WORK( N2+1 ), SVA( N2+1 ), MVL, |
+ 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 ) |
+ WORK( N+1 ), LWORK-N, IERR ) |
* |
* |
CALL DGSVJ0( JOBV, M-N4, N2-N4, A( N4+1, N4+1 ), LDA, |
CALL DGSVJ0( JOBV, M-N4, N2-N4, A( N4+1, N4+1 ), LDA, |
+ WORK( N4+1 ), SVA( N4+1 ), MVL, |
+ 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 ) |
+ WORK( N+1 ), LWORK-N, IERR ) |
* |
* |
CALL DGSVJ0( JOBV, M, N4, A, LDA, WORK, SVA, MVL, V, LDV, |
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 ) |
+ IERR ) |
* |
* |
CALL DGSVJ1( JOBV, M, N2, N4, A, LDA, WORK, SVA, MVL, V, |
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 ) |
+ LWORK-N, IERR ) |
* |
* |
* |
* |
Line 670
|
Line 670
|
* |
* |
* |
* |
CALL DGSVJ0( JOBV, N4, N4, A, LDA, WORK, SVA, MVL, V, LDV, |
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 ) |
+ IERR ) |
* |
* |
CALL DGSVJ0( JOBV, N2, N4, A( 1, N4+1 ), LDA, WORK( N4+1 ), |
CALL DGSVJ0( JOBV, N2, N4, A( 1, N4+1 ), LDA, WORK( N4+1 ), |
+ SVA( N4+1 ), MVL, V( N4*q+1, N4+1 ), LDV, |
+ 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 ) |
+ IERR ) |
* |
* |
CALL DGSVJ1( JOBV, N2, N2, N4, A, LDA, WORK, SVA, MVL, V, |
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 ) |
+ LWORK-N, IERR ) |
* |
* |
CALL DGSVJ0( JOBV, N2+N4, N4, A( 1, N2+1 ), LDA, |
CALL DGSVJ0( JOBV, N2+N4, N4, A( 1, N2+1 ), LDA, |
+ WORK( N2+1 ), SVA( N2+1 ), MVL, |
+ 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 ) |
+ WORK( N+1 ), LWORK-N, IERR ) |
|
|
END IF |
END IF |
Line 753
|
Line 753
|
SVA( p ) = DNRM2( M, A( 1, p ), 1 )*WORK( p ) |
SVA( p ) = DNRM2( M, A( 1, p ), 1 )*WORK( p ) |
ELSE |
ELSE |
TEMP1 = ZERO |
TEMP1 = ZERO |
AAPP = ZERO |
AAPP = ONE |
CALL DLASSQ( M, A( 1, p ), 1, TEMP1, AAPP ) |
CALL DLASSQ( M, A( 1, p ), 1, TEMP1, AAPP ) |
SVA( p ) = TEMP1*DSQRT( AAPP )*WORK( p ) |
SVA( p ) = TEMP1*DSQRT( AAPP )*WORK( p ) |
END IF |
END IF |
Line 841
|
Line 841
|
+ FASTR ) |
+ FASTR ) |
SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, |
SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, |
+ ONE+T*APOAQ*AAPQ ) ) |
+ ONE+T*APOAQ*AAPQ ) ) |
AAPP = AAPP*DSQRT( ONE-T*AQOAP* |
AAPP = AAPP*DSQRT( DMAX1( ZERO, |
+ AAPQ ) |
+ ONE-T*AQOAP*AAPQ ) ) |
MXSINJ = DMAX1( MXSINJ, DABS( T ) ) |
MXSINJ = DMAX1( MXSINJ, DABS( T ) ) |
* |
* |
ELSE |
ELSE |
Line 989
|
Line 989
|
+ WORK( q ) |
+ WORK( q ) |
ELSE |
ELSE |
T = ZERO |
T = ZERO |
AAQQ = ZERO |
AAQQ = ONE |
CALL DLASSQ( M, A( 1, q ), 1, T, |
CALL DLASSQ( M, A( 1, q ), 1, T, |
+ AAQQ ) |
+ AAQQ ) |
SVA( q ) = T*DSQRT( AAQQ )*WORK( q ) |
SVA( q ) = T*DSQRT( AAQQ )*WORK( q ) |
Line 1002
|
Line 1002
|
+ WORK( p ) |
+ WORK( p ) |
ELSE |
ELSE |
T = ZERO |
T = ZERO |
AAPP = ZERO |
AAPP = ONE |
CALL DLASSQ( M, A( 1, p ), 1, T, |
CALL DLASSQ( M, A( 1, p ), 1, T, |
+ AAPP ) |
+ AAPP ) |
AAPP = T*DSQRT( AAPP )*WORK( p ) |
AAPP = T*DSQRT( AAPP )*WORK( p ) |
Line 1164
|
Line 1164
|
MXSINJ = DMAX1( MXSINJ, DABS( SN ) ) |
MXSINJ = DMAX1( MXSINJ, DABS( SN ) ) |
SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, |
SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, |
+ ONE+T*APOAQ*AAPQ ) ) |
+ ONE+T*APOAQ*AAPQ ) ) |
AAPP = AAPP*DSQRT( ONE-T*AQOAP* |
AAPP = AAPP*DSQRT( DMAX1( ZERO, |
+ AAPQ ) |
+ ONE-T*AQOAP*AAPQ ) ) |
* |
* |
APOAQ = WORK( p ) / WORK( q ) |
APOAQ = WORK( p ) / WORK( q ) |
AQOAP = WORK( q ) / WORK( p ) |
AQOAP = WORK( q ) / WORK( p ) |
Line 1316
|
Line 1316
|
+ WORK( q ) |
+ WORK( q ) |
ELSE |
ELSE |
T = ZERO |
T = ZERO |
AAQQ = ZERO |
AAQQ = ONE |
CALL DLASSQ( M, A( 1, q ), 1, T, |
CALL DLASSQ( M, A( 1, q ), 1, T, |
+ AAQQ ) |
+ AAQQ ) |
SVA( q ) = T*DSQRT( AAQQ )*WORK( q ) |
SVA( q ) = T*DSQRT( AAQQ )*WORK( q ) |
Line 1329
|
Line 1329
|
+ WORK( p ) |
+ WORK( p ) |
ELSE |
ELSE |
T = ZERO |
T = ZERO |
AAPP = ZERO |
AAPP = ONE |
CALL DLASSQ( M, A( 1, p ), 1, T, |
CALL DLASSQ( M, A( 1, p ), 1, T, |
+ AAPP ) |
+ AAPP ) |
AAPP = T*DSQRT( AAPP )*WORK( p ) |
AAPP = T*DSQRT( AAPP )*WORK( p ) |
Line 1395
|
Line 1395
|
SVA( N ) = DNRM2( M, A( 1, N ), 1 )*WORK( N ) |
SVA( N ) = DNRM2( M, A( 1, N ), 1 )*WORK( N ) |
ELSE |
ELSE |
T = ZERO |
T = ZERO |
AAPP = ZERO |
AAPP = ONE |
CALL DLASSQ( M, A( 1, N ), 1, T, AAPP ) |
CALL DLASSQ( M, A( 1, N ), 1, T, AAPP ) |
SVA( N ) = T*DSQRT( AAPP )*WORK( N ) |
SVA( N ) = T*DSQRT( AAPP )*WORK( N ) |
END IF |
END IF |
Line 1446
|
Line 1446
|
END IF |
END IF |
IF( SVA( p ).NE.ZERO ) THEN |
IF( SVA( p ).NE.ZERO ) THEN |
N4 = N4 + 1 |
N4 = N4 + 1 |
IF( SVA( p )*SCALE.GT.SFMIN )N2 = N2 + 1 |
IF( SVA( p )*SKL.GT.SFMIN )N2 = N2 + 1 |
END IF |
END IF |
5991 CONTINUE |
5991 CONTINUE |
IF( SVA( N ).NE.ZERO ) THEN |
IF( SVA( N ).NE.ZERO ) THEN |
N4 = N4 + 1 |
N4 = N4 + 1 |
IF( SVA( N )*SCALE.GT.SFMIN )N2 = N2 + 1 |
IF( SVA( N )*SKL.GT.SFMIN )N2 = N2 + 1 |
END IF |
END IF |
* |
* |
* Normalize the left singular vectors. |
* Normalize the left singular vectors. |
Line 1478
|
Line 1478
|
END IF |
END IF |
* |
* |
* Undo scaling, if necessary (and possible). |
* Undo scaling, if necessary (and possible). |
IF( ( ( SCALE.GT.ONE ) .AND. ( SVA( 1 ).LT.( BIG / |
IF( ( ( SKL.GT.ONE ) .AND. ( SVA( 1 ).LT.( BIG / |
+ SCALE ) ) ) .OR. ( ( SCALE.LT.ONE ) .AND. ( SVA( N2 ).GT. |
+ SKL) ) ) .OR. ( ( SKL.LT.ONE ) .AND. ( SVA( N2 ).GT. |
+ ( SFMIN / SCALE ) ) ) ) THEN |
+ ( SFMIN / SKL) ) ) ) THEN |
DO 2400 p = 1, N |
DO 2400 p = 1, N |
SVA( p ) = SCALE*SVA( p ) |
SVA( p ) = SKL*SVA( p ) |
2400 CONTINUE |
2400 CONTINUE |
SCALE = ONE |
SKL= ONE |
END IF |
END IF |
* |
* |
WORK( 1 ) = SCALE |
WORK( 1 ) = SKL |
* The singular values of A are SCALE*SVA(1:N). If SCALE.NE.ONE |
* 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 |
* then some of the singular values may overflow or underflow and |
* the spectrum is given in this factored representation. |
* the spectrum is given in this factored representation. |
* |
* |