diff options
-rw-r--r-- | SRC/dgesvj.f | 14 | ||||
-rw-r--r-- | SRC/dgsvj0.f | 14 | ||||
-rw-r--r-- | SRC/dgsvj1.f | 7 | ||||
-rw-r--r-- | SRC/sgejsv.f | 79 | ||||
-rw-r--r-- | SRC/sgesvj.f | 16 | ||||
-rw-r--r-- | SRC/sgsvj0.f | 5 |
6 files changed, 81 insertions, 54 deletions
diff --git a/SRC/dgesvj.f b/SRC/dgesvj.f index d48975e3..a8e154fc 100644 --- a/SRC/dgesvj.f +++ b/SRC/dgesvj.f @@ -1,11 +1,11 @@ SUBROUTINE DGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, + LDV, WORK, LWORK, INFO ) * -* -- LAPACK routine (version 3.3.0) -- +* -- LAPACK routine (version 3.3.1) -- * * -- Contributed by Zlatko Drmac of the University of Zagreb and -- * -- Kresimir Veselic of the Fernuniversitaet Hagen -- -* November 2010 +* January 2011 * * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- @@ -133,7 +133,7 @@ * referenced * * M (input) INTEGER -* The number of rows of the input matrix A. M >= 0. +* The number of rows of the input matrix A. 1/DLAMCH('E') > M >= 0. * * N (input) INTEGER * The number of columns of the input matrix A. @@ -383,7 +383,7 @@ ROOTTOL = DSQRT( TOL ) * IF( DBLE( M )*EPSLN.GE.ONE ) THEN - INFO = -5 + INFO = -4 CALL XERBLA( 'DGESVJ', -INFO ) RETURN END IF @@ -824,8 +824,7 @@ * AQOAP = AAQQ / AAPP APOAQ = AAPP / AAQQ - THETA = -HALF*DABS( AQOAP-APOAQ ) / - + AAPQ + THETA = -HALF*DABS(AQOAP-APOAQ)/AAPQ * IF( DABS( THETA ).GT.BIGTHETA ) THEN * @@ -1131,8 +1130,7 @@ * AQOAP = AAQQ / AAPP APOAQ = AAPP / AAQQ - THETA = -HALF*DABS( AQOAP-APOAQ ) / - + AAPQ + THETA = -HALF*DABS(AQOAP-APOAQ)/AAPQ IF( AAQQ.GT.AAPP0 )THETA = -THETA * IF( DABS( THETA ).GT.BIGTHETA ) THEN diff --git a/SRC/dgsvj0.f b/SRC/dgsvj0.f index e03aeb63..c26f929f 100644 --- a/SRC/dgsvj0.f +++ b/SRC/dgsvj0.f @@ -1,11 +1,11 @@ SUBROUTINE DGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, + SFMIN, TOL, NSWEEP, WORK, LWORK, INFO ) * -* -- LAPACK routine (version 3.3.0) -- +* -- LAPACK routine (version 3.3.1) -- * * -- Contributed by Zlatko Drmac of the University of Zagreb and -- * -- Kresimir Veselic of the Fernuniversitaet Hagen -- -* November 2010 +* January 2011 * * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- @@ -16,6 +16,7 @@ * eigenvalue problems Hx = lambda M x, H M x = lambda x with H, M > 0. * IMPLICIT NONE +* .. * .. Scalar Arguments .. INTEGER INFO, LDA, LDV, LWORK, M, MV, N, NSWEEP DOUBLE PRECISION EPS, SFMIN, TOL @@ -173,6 +174,8 @@ * .. * .. Executable Statements .. * +* Test the input parameters. +* APPLV = LSAME( JOBV, 'A' ) RSVEC = LSAME( JOBV, 'V' ) IF( .NOT.( RSVEC .OR. APPLV .OR. LSAME( JOBV, 'N' ) ) ) THEN @@ -219,7 +222,6 @@ BIGTHETA = ONE / ROOTEPS ROOTTOL = DSQRT( TOL ) * -* * -#- Row-cyclic Jacobi SVD algorithm with column pivoting -#- * EMPTSW = ( N*( N-1 ) ) / 2 @@ -376,8 +378,7 @@ * AQOAP = AAQQ / AAPP APOAQ = AAPP / AAQQ - THETA = -HALF*DABS( AQOAP-APOAQ ) / - + AAPQ + THETA = -HALF*DABS( AQOAP-APOAQ )/AAPQ * IF( DABS( THETA ).GT.BIGTHETA ) THEN * @@ -676,8 +677,7 @@ * AQOAP = AAQQ / AAPP APOAQ = AAPP / AAQQ - THETA = -HALF*DABS( AQOAP-APOAQ ) / - + AAPQ + THETA = -HALF*DABS( AQOAP-APOAQ )/AAPQ IF( AAQQ.GT.AAPP0 )THETA = -THETA * IF( DABS( THETA ).GT.BIGTHETA ) THEN diff --git a/SRC/dgsvj1.f b/SRC/dgsvj1.f index 1a51f15c..89704e73 100644 --- a/SRC/dgsvj1.f +++ b/SRC/dgsvj1.f @@ -1,11 +1,11 @@ SUBROUTINE DGSVJ1( JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV, + EPS, SFMIN, TOL, NSWEEP, WORK, LWORK, INFO ) * -* -- LAPACK routine (version 3.3.0) -- +* -- LAPACK routine (version 3.3.1) -- * * -- Contributed by Zlatko Drmac of the University of Zagreb and -- * -- Kresimir Veselic of the Fernuniversitaet Hagen -- -* November 2010 +* January 2011 * * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- @@ -375,8 +375,7 @@ * AQOAP = AAQQ / AAPP APOAQ = AAPP / AAQQ - THETA = -HALF*DABS( AQOAP-APOAQ ) / - + AAPQ + THETA = -HALF*DABS(AQOAP-APOAQ) / AAPQ IF( AAQQ.GT.AAPP0 )THETA = -THETA IF( DABS( THETA ).GT.BIGTHETA ) THEN diff --git a/SRC/sgejsv.f b/SRC/sgejsv.f index e2138cff..aaef93f4 100644 --- a/SRC/sgejsv.f +++ b/SRC/sgejsv.f @@ -2,11 +2,11 @@ & M, N, A, LDA, SVA, U, LDU, V, LDV, & WORK, LWORK, IWORK, INFO ) * -* -- LAPACK routine (version 3.3.0) -- +* -- LAPACK routine (version 3.3.1) -- * * -- Contributed by Zlatko Drmac of the University of Zagreb and -- * -- Kresimir Veselic of the Fernuniversitaet Hagen -- -* November 2010 +* January 2011 * * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- @@ -31,6 +31,7 @@ * * Purpose * ======= +* * SGEJSV computes the singular value decomposition (SVD) of a real M-by-N * matrix [A], where M >= N. The SVD of [A] is written as * @@ -253,33 +254,50 @@ * LWORK depends on the job: * * If only SIGMA is needed ( JOBU.EQ.'N', JOBV.EQ.'N' ) and -* -> .. no scaled condition estimate required ( JOBE.EQ.'N'): +* -> .. no scaled condition estimate required (JOBE.EQ.'N'): * LWORK >= max(2*M+N,4*N+1,7). This is the minimal requirement. -* For optimal performance (blocked code) the optimal value +* ->> For optimal performance (blocked code) the optimal value * is LWORK >= max(2*M+N,3*N+(N+1)*NB,7). Here NB is the optimal -* block size for xGEQP3/xGEQRF. +* block size for DGEQP3 and DGEQRF. +* In general, optimal LWORK is computed as +* LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DGEQRF), 7). * -> .. an estimate of the scaled condition number of A is * required (JOBA='E', 'G'). In this case, LWORK is the maximum -* of the above and N*N+4*N, i.e. LWORK >= max(2*M+N,N*N+4N,7). +* of the above and N*N+4*N, i.e. LWORK >= max(2*M+N,N*N+4*N,7). +* ->> For optimal performance (blocked code) the optimal value +* is LWORK >= max(2*M+N,3*N+(N+1)*NB, N*N+4*N, 7). +* In general, the optimal length LWORK is computed as +* LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DGEQRF), +* N+N*N+LWORK(DPOCON),7). * * If SIGMA and the right singular vectors are needed (JOBV.EQ.'V'), -* -> the minimal requirement is LWORK >= max(2*N+M,7). -* -> For optimal performance, LWORK >= max(2*N+M,2*N+N*NB,7), -* where NB is the optimal block size. +* -> the minimal requirement is LWORK >= max(2*M+N,4*N+1,7). +* -> For optimal performance, LWORK >= max(2*M+N,3*N+(N+1)*NB,7), +* where NB is the optimal block size for DGEQP3, DGEQRF, DGELQ, +* DORMLQ. In general, the optimal length LWORK is computed as +* LWORK >= max(2*M+N,N+LWORK(DGEQP3), N+LWORK(DPOCON), +* N+LWORK(DGELQ), 2*N+LWORK(DGEQRF), N+LWORK(DORMLQ)). * * If SIGMA and the left singular vectors are needed -* -> the minimal requirement is LWORK >= max(2*N+M,7). -* -> For optimal performance, LWORK >= max(2*N+M,2*N+N*NB,7), -* where NB is the optimal block size. -* -* If full SVD is needed ( JOBU.EQ.'U' or 'F', JOBV.EQ.'V' ) and -* -> .. the singular vectors are computed without explicit -* accumulation of the Jacobi rotations, LWORK >= 6*N+2*N*N -* -> .. in the iterative part, the Jacobi rotations are -* explicitly accumulated (option, see the description of JOBV), -* then the minimal requirement is LWORK >= max(M+3*N+N*N,7). -* For better performance, if NB is the optimal block size, -* LWORK >= max(3*N+N*N+M,3*N+N*N+N*NB,7). +* -> the minimal requirement is LWORK >= max(2*M+N,4*N+1,7). +* -> For optimal performance: +* if JOBU.EQ.'U' :: LWORK >= max(2*M+N,3*N+(N+1)*NB,7), +* if JOBU.EQ.'F' :: LWORK >= max(2*M+N,3*N+(N+1)*NB,N+M*NB,7), +* where NB is the optimal block size for DGEQP3, DGEQRF, DORMQR. +* In general, the optimal length LWORK is computed as +* LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DPOCON), +* 2*N+LWORK(DGEQRF), N+LWORK(DORMQR)). +* Here LWORK(DORMQR) equals N*NB (for JOBU.EQ.'U') or +* M*NB (for JOBU.EQ.'F'). +* +* If the full SVD is needed: (JOBU.EQ.'U' or JOBU.EQ.'F') and +* -> if JOBV.EQ.'V' +* the minimal requirement is LWORK >= max(2*M+N,6*N+2*N*N). +* -> if JOBV.EQ.'J' the minimal requirement is +* LWORK >= max(2*M+N, 4*N+N*N,2*N+N*N+6). +* -> For optimal performance, LWORK should be additionally +* larger than N+M*NB, where NB is the optimal block size +* for DORMQR. * * IWORK (workspace/output) INTEGER array, dimension M+3*N. * On exit, @@ -439,12 +457,16 @@ INFO = - 14 ELSE IF ( (.NOT.(LSVEC .OR. RSVEC .OR. ERREST).AND. & (LWORK .LT. MAX0(7,4*N+1,2*M+N))) .OR. - & (.NOT.(LSVEC .OR. LSVEC) .AND. ERREST .AND. + & (.NOT.(LSVEC .OR. RSVEC) .AND. ERREST .AND. & (LWORK .LT. MAX0(7,4*N+N*N,2*M+N))) .OR. - & (LSVEC .AND. (.NOT.RSVEC) .AND. (LWORK .LT. MAX0(7,2*N+M))) .OR. - & (RSVEC .AND. (.NOT.LSVEC) .AND. (LWORK .LT. MAX0(7,2*N+M))) .OR. - & (LSVEC .AND. RSVEC .AND. .NOT.JRACC .AND. (LWORK.LT.6*N+2*N*N)) - & .OR. (LSVEC.AND.RSVEC.AND.JRACC.AND.LWORK.LT.MAX0(7,M+3*N+N*N))) + & (LSVEC .AND. (.NOT.RSVEC) .AND. (LWORK .LT. MAX0(7,2*M+N,4*N+1))) + & .OR. + & (RSVEC .AND. (.NOT.LSVEC) .AND. (LWORK .LT. MAX0(7,2*M+N,4*N+1))) + & .OR. + & (LSVEC .AND. RSVEC .AND. (.NOT.JRACC) .AND. + & (LWORK.LT.MAX0(2*M+N,6*N+2*N*N))) + & .OR. (LSVEC .AND. RSVEC .AND. JRACC .AND. + & LWORK.LT.MAX0(2*M+N,4*N+N*N,2*N+N*N+6))) & THEN INFO = - 17 ELSE @@ -455,6 +477,7 @@ IF ( INFO .NE. 0 ) THEN * #:( CALL XERBLA( 'SGEJSV', - INFO ) + RETURN END IF * * Quick return for void matrix (Y3K safe) @@ -476,6 +499,7 @@ SFMIN = SLAMCH('SafeMinimum') SMALL = SFMIN / EPSLN BIG = SLAMCH('O') +* BIG = ONE / SFMIN * * Initialize SVA(1:N) = diag( ||A e_i||_2 )_1^N * @@ -535,6 +559,7 @@ END IF IWORK(1) = 0 IWORK(2) = 0 + IWORK(3) = 0 RETURN END IF * @@ -1028,7 +1053,7 @@ 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, INFO ) + & LDU, WORK(N+1), LWORK-N, INFO ) SCALEM = WORK(N+1) NUMRANK = NINT(WORK(N+2)) IF ( NR .LT. N ) THEN diff --git a/SRC/sgesvj.f b/SRC/sgesvj.f index 01a3d248..f5fa3978 100644 --- a/SRC/sgesvj.f +++ b/SRC/sgesvj.f @@ -1,11 +1,11 @@ SUBROUTINE SGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, + LDV, WORK, LWORK, INFO ) * -* -- LAPACK routine (version 3.3.0) -- +* -- LAPACK routine (version 3.3.1) -- * * -- Contributed by Zlatko Drmac of the University of Zagreb and -- * -- Kresimir Veselic of the Fernuniversitaet Hagen -- -* November 2010 +* January 2011 * * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- @@ -133,7 +133,7 @@ * referenced * * M (input) INTEGER -* The number of rows of the input matrix A. M >= 0. +* The number of rows of the input matrix A. 1/SLAMCH('E') > M >= 0. * * N (input) INTEGER * The number of columns of the input matrix A. @@ -241,7 +241,8 @@ * Jacobi rotation angles in the last sweep. It can be * useful for a post festum analysis. * -* LWORK length of WORK, WORK >= MAX(6,M+N) +* LWORK (input) INTEGER +* length of WORK, WORK >= MAX(6,M+N) * * INFO (output) INTEGER * = 0 : successful exit. @@ -249,6 +250,7 @@ * > 0 : SGESVJ did not converge in the maximal allowed number (30) * of sweeps. The output may still be useful. See the * description of WORK. +* * ===================================================================== * * .. Local Parameters .. @@ -278,6 +280,7 @@ INTRINSIC ABS, AMAX1, AMIN1, FLOAT, MIN0, SIGN, SQRT * .. * .. External Functions .. +* .. * from BLAS REAL SDOT, SNRM2 EXTERNAL SDOT, SNRM2 @@ -290,6 +293,7 @@ EXTERNAL LSAME * .. * .. External Subroutines .. +* .. * from BLAS EXTERNAL SAXPY, SCOPY, SROTM, SSCAL, SSWAP * from LAPACK @@ -370,6 +374,7 @@ ROOTSFMIN = SQRT( SFMIN ) SMALL = SFMIN / EPSLN BIG = SLAMCH( 'Overflow' ) +* BIG = ONE / SFMIN ROOTBIG = ONE / ROOTSFMIN LARGE = BIG / SQRT( FLOAT( M*N ) ) BIGTHETA = ONE / ROOTEPS @@ -378,7 +383,7 @@ ROOTTOL = SQRT( TOL ) * IF( FLOAT( M )*EPSLN.GE.ONE ) THEN - INFO = -5 + INFO = -4 CALL XERBLA( 'SGESVJ', -INFO ) RETURN END IF @@ -689,6 +694,7 @@ * .. Row-cyclic pivot strategy with de Rijk's pivoting .. * DO 1993 i = 1, NSWEEP +* * .. go go go ... * MXAAPQ = ZERO diff --git a/SRC/sgsvj0.f b/SRC/sgsvj0.f index 1279162d..07c4d030 100644 --- a/SRC/sgsvj0.f +++ b/SRC/sgsvj0.f @@ -1,11 +1,11 @@ SUBROUTINE SGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, + SFMIN, TOL, NSWEEP, WORK, LWORK, INFO ) * -* -- LAPACK routine (version 3.3.0) -- +* -- LAPACK routine (version 3.3.1) -- * * -- Contributed by Zlatko Drmac of the University of Zagreb and -- * -- Kresimir Veselic of the Fernuniversitaet Hagen -- -* November 2010 +* January 2011 * * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- @@ -222,7 +222,6 @@ BIGTHETA = ONE / ROOTEPS ROOTTOL = SQRT( TOL ) * -* * .. Row-cyclic Jacobi SVD algorithm with column pivoting .. * EMPTSW = ( N*( N-1 ) ) / 2 |