diff --git a/lapack-netlib/SRC/cbdsqr.f b/lapack-netlib/SRC/cbdsqr.f index cf1459ad22..fd6d4731b7 100644 --- a/lapack-netlib/SRC/cbdsqr.f +++ b/lapack-netlib/SRC/cbdsqr.f @@ -5,7 +5,6 @@ * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> \htmlonly *> Download CBDSQR + dependencies *> *> [TGZ] @@ -13,7 +12,6 @@ *> [ZIP] *> *> [TXT] -*> \endhtmlonly * * Definition: * =========== @@ -166,7 +164,9 @@ *> *> \param[out] RWORK *> \verbatim -*> RWORK is REAL array, dimension (4*N) +*> RWORK is REAL array, dimension (LRWORK) +*> LRWORK = 4*N, if NCVT = NRU = NCC = 0, and +*> LRWORK = 4*(N-1), otherwise *> \endverbatim *> *> \param[out] INFO @@ -230,6 +230,7 @@ * ===================================================================== SUBROUTINE CBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, $ LDU, C, LDC, RWORK, INFO ) + IMPLICIT NONE * * -- LAPACK computational routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -279,7 +280,8 @@ SUBROUTINE CBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, EXTERNAL LSAME, SLAMCH * .. * .. External Subroutines .. - EXTERNAL CLASR, CSROT, CSSCAL, CSWAP, SLARTG, SLAS2, + EXTERNAL CLASR, CSROT, CSSCAL, CSWAP, SLARTG, + $ SLAS2, $ SLASQ1, SLASV2, XERBLA * .. * .. Intrinsic Functions .. @@ -328,9 +330,10 @@ SUBROUTINE CBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, IF( .NOT.ROTATE ) THEN CALL SLASQ1( N, D, E, RWORK, INFO ) * -* If INFO equals 2, dqds didn't finish, try to finish +* If the dqds algorithm failed, try to finish with +* the standard QR algorithm * - IF( INFO .NE. 2 ) RETURN + IF( INFO .EQ. 0 ) RETURN INFO = 0 END IF * @@ -360,10 +363,12 @@ SUBROUTINE CBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, * Update singular vectors if desired * IF( NRU.GT.0 ) - $ CALL CLASR( 'R', 'V', 'F', NRU, N, RWORK( 1 ), RWORK( N ), + $ CALL CLASR( 'R', 'V', 'F', NRU, N, RWORK( 1 ), + $ RWORK( N ), $ U, LDU ) IF( NCC.GT.0 ) - $ CALL CLASR( 'L', 'V', 'F', N, NCC, RWORK( 1 ), RWORK( N ), + $ CALL CLASR( 'L', 'V', 'F', N, NCC, RWORK( 1 ), + $ RWORK( N ), $ C, LDC ) END IF * @@ -400,12 +405,14 @@ SUBROUTINE CBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, 40 CONTINUE 50 CONTINUE SMINOA = SMINOA / SQRT( REAL( N ) ) - THRESH = MAX( TOL*SMINOA, MAXITR*(N*(N*UNFL)) ) + THRESH = MAX( TOL*SMINOA, + $ REAL(MAXITR)*(REAL(N)*(REAL(N)*UNFL)) ) ELSE * * Absolute accuracy desired * - THRESH = MAX( ABS( TOL )*SMAX, MAXITR*(N*(N*UNFL)) ) + THRESH = MAX( ABS( TOL )*SMAX, + $ REAL(MAXITR)*(REAL(N)*(REAL(N)*UNFL)) ) END IF * * Prepare for main iteration loop for the singular values @@ -487,7 +494,8 @@ SUBROUTINE CBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, $ CALL CSROT( NCVT, VT( M-1, 1 ), LDVT, VT( M, 1 ), LDVT, $ COSR, SINR ) IF( NRU.GT.0 ) - $ CALL CSROT( NRU, U( 1, M-1 ), 1, U( 1, M ), 1, COSL, SINL ) + $ CALL CSROT( NRU, U( 1, M-1 ), 1, U( 1, M ), 1, COSL, + $ SINL ) IF( NCC.GT.0 ) $ CALL CSROT( NCC, C( M-1, 1 ), LDC, C( M, 1 ), LDC, COSL, $ SINL ) @@ -576,7 +584,7 @@ SUBROUTINE CBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, * Compute shift. First, test if shifting would ruin relative * accuracy, and if so set the shift to zero. * - IF( TOL.GE.ZERO .AND. N*TOL*( SMIN / SMAX ).LE. + IF( TOL.GE.ZERO .AND. REAL( N )*TOL*( SMIN / SMAX ).LE. $ MAX( EPS, HNDRTH*TOL ) ) THEN * * Use a zero shift to avoid loss of relative accuracy @@ -620,7 +628,8 @@ SUBROUTINE CBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, CALL SLARTG( D( I )*CS, E( I ), CS, SN, R ) IF( I.GT.LL ) $ E( I-1 ) = OLDSN*R - CALL SLARTG( OLDCS*R, D( I+1 )*SN, OLDCS, OLDSN, D( I ) ) + CALL SLARTG( OLDCS*R, D( I+1 )*SN, OLDCS, OLDSN, + $ D( I ) ) RWORK( I-LL+1 ) = CS RWORK( I-LL+1+NM1 ) = SN RWORK( I-LL+1+NM12 ) = OLDCS @@ -636,10 +645,12 @@ SUBROUTINE CBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, $ CALL CLASR( 'L', 'V', 'F', M-LL+1, NCVT, RWORK( 1 ), $ RWORK( N ), VT( LL, 1 ), LDVT ) IF( NRU.GT.0 ) - $ CALL CLASR( 'R', 'V', 'F', NRU, M-LL+1, RWORK( NM12+1 ), + $ CALL CLASR( 'R', 'V', 'F', NRU, M-LL+1, + $ RWORK( NM12+1 ), $ RWORK( NM13+1 ), U( 1, LL ), LDU ) IF( NCC.GT.0 ) - $ CALL CLASR( 'L', 'V', 'F', M-LL+1, NCC, RWORK( NM12+1 ), + $ CALL CLASR( 'L', 'V', 'F', M-LL+1, NCC, + $ RWORK( NM12+1 ), $ RWORK( NM13+1 ), C( LL, 1 ), LDC ) * * Test convergence @@ -658,7 +669,8 @@ SUBROUTINE CBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, CALL SLARTG( D( I )*CS, E( I-1 ), CS, SN, R ) IF( I.LT.M ) $ E( I ) = OLDSN*R - CALL SLARTG( OLDCS*R, D( I-1 )*SN, OLDCS, OLDSN, D( I ) ) + CALL SLARTG( OLDCS*R, D( I-1 )*SN, OLDCS, OLDSN, + $ D( I ) ) RWORK( I-LL ) = CS RWORK( I-LL+NM1 ) = -SN RWORK( I-LL+NM12 ) = OLDCS @@ -671,7 +683,8 @@ SUBROUTINE CBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, * Update singular vectors * IF( NCVT.GT.0 ) - $ CALL CLASR( 'L', 'V', 'B', M-LL+1, NCVT, RWORK( NM12+1 ), + $ CALL CLASR( 'L', 'V', 'B', M-LL+1, NCVT, + $ RWORK( NM12+1 ), $ RWORK( NM13+1 ), VT( LL, 1 ), LDVT ) IF( NRU.GT.0 ) $ CALL CLASR( 'R', 'V', 'B', NRU, M-LL+1, RWORK( 1 ), @@ -726,10 +739,12 @@ SUBROUTINE CBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, $ CALL CLASR( 'L', 'V', 'F', M-LL+1, NCVT, RWORK( 1 ), $ RWORK( N ), VT( LL, 1 ), LDVT ) IF( NRU.GT.0 ) - $ CALL CLASR( 'R', 'V', 'F', NRU, M-LL+1, RWORK( NM12+1 ), + $ CALL CLASR( 'R', 'V', 'F', NRU, M-LL+1, + $ RWORK( NM12+1 ), $ RWORK( NM13+1 ), U( 1, LL ), LDU ) IF( NCC.GT.0 ) - $ CALL CLASR( 'L', 'V', 'F', M-LL+1, NCC, RWORK( NM12+1 ), + $ CALL CLASR( 'L', 'V', 'F', M-LL+1, NCC, + $ RWORK( NM12+1 ), $ RWORK( NM13+1 ), C( LL, 1 ), LDC ) * * Test convergence @@ -776,7 +791,8 @@ SUBROUTINE CBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, * Update singular vectors if desired * IF( NCVT.GT.0 ) - $ CALL CLASR( 'L', 'V', 'B', M-LL+1, NCVT, RWORK( NM12+1 ), + $ CALL CLASR( 'L', 'V', 'B', M-LL+1, NCVT, + $ RWORK( NM12+1 ), $ RWORK( NM13+1 ), VT( LL, 1 ), LDVT ) IF( NRU.GT.0 ) $ CALL CLASR( 'R', 'V', 'B', NRU, M-LL+1, RWORK( 1 ), @@ -795,6 +811,12 @@ SUBROUTINE CBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, * 160 CONTINUE DO 170 I = 1, N + IF( D( I ).EQ.ZERO ) THEN +* +* Avoid -ZERO +* + D( I ) = ZERO + END IF IF( D( I ).LT.ZERO ) THEN D( I ) = -D( I ) * @@ -832,7 +854,8 @@ SUBROUTINE CBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, IF( NRU.GT.0 ) $ CALL CSWAP( NRU, U( 1, ISUB ), 1, U( 1, N+1-I ), 1 ) IF( NCC.GT.0 ) - $ CALL CSWAP( NCC, C( ISUB, 1 ), LDC, C( N+1-I, 1 ), LDC ) + $ CALL CSWAP( NCC, C( ISUB, 1 ), LDC, C( N+1-I, 1 ), + $ LDC ) END IF 190 CONTINUE GO TO 220 diff --git a/lapack-netlib/SRC/dbdsqr.f b/lapack-netlib/SRC/dbdsqr.f index bc697a0074..99a12aa79a 100644 --- a/lapack-netlib/SRC/dbdsqr.f +++ b/lapack-netlib/SRC/dbdsqr.f @@ -5,7 +5,6 @@ * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> \htmlonly *> Download DBDSQR + dependencies *> *> [TGZ] @@ -13,7 +12,6 @@ *> [ZIP] *> *> [TXT] -*> \endhtmlonly * * Definition: * =========== @@ -166,7 +164,9 @@ *> *> \param[out] WORK *> \verbatim -*> WORK is DOUBLE PRECISION array, dimension (4*(N-1)) +*> WORK is DOUBLE PRECISION array, dimension (LWORK) +*> LWORK = 4*N, if NCVT = NRU = NCC = 0, and +*> LWORK = 4*(N-1), otherwise *> \endverbatim *> *> \param[out] INFO @@ -181,7 +181,7 @@ *> iterations (in inner while loop) *> = 3, termination criterion of outer while loop not met *> (program created more than N unreduced blocks) -*> else NCVT = NRU = NCC = 0, +*> else, *> the algorithm did not converge; D and E contain the *> elements of a bidiagonal matrix which is orthogonally *> similar to the input matrix B; if INFO = i, i @@ -233,11 +233,12 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup auxOTHERcomputational +*> \ingroup bdsqr * * ===================================================================== SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, $ LDU, C, LDC, WORK, INFO ) + IMPLICIT NONE * * -- LAPACK computational routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -287,7 +288,8 @@ SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, EXTERNAL LSAME, DLAMCH * .. * .. External Subroutines .. - EXTERNAL DLARTG, DLAS2, DLASQ1, DLASR, DLASV2, DROT, + EXTERNAL DLARTG, DLAS2, DLASQ1, DLASR, DLASV2, + $ DROT, $ DSCAL, DSWAP, XERBLA * .. * .. Intrinsic Functions .. @@ -336,9 +338,10 @@ SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, IF( .NOT.ROTATE ) THEN CALL DLASQ1( N, D, E, WORK, INFO ) * -* If INFO equals 2, dqds didn't finish, try to finish +* If the dqds algorithm failed, try to finish with +* the standard QR algorithm * - IF( INFO .NE. 2 ) RETURN + IF( INFO .EQ. 0 ) RETURN INFO = 0 END IF * @@ -368,10 +371,12 @@ SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, * Update singular vectors if desired * IF( NRU.GT.0 ) - $ CALL DLASR( 'R', 'V', 'F', NRU, N, WORK( 1 ), WORK( N ), U, + $ CALL DLASR( 'R', 'V', 'F', NRU, N, WORK( 1 ), WORK( N ), + $ U, $ LDU ) IF( NCC.GT.0 ) - $ CALL DLASR( 'L', 'V', 'F', N, NCC, WORK( 1 ), WORK( N ), C, + $ CALL DLASR( 'L', 'V', 'F', N, NCC, WORK( 1 ), WORK( N ), + $ C, $ LDC ) END IF * @@ -493,10 +498,12 @@ SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, * Compute singular vectors, if desired * IF( NCVT.GT.0 ) - $ CALL DROT( NCVT, VT( M-1, 1 ), LDVT, VT( M, 1 ), LDVT, COSR, + $ CALL DROT( NCVT, VT( M-1, 1 ), LDVT, VT( M, 1 ), LDVT, + $ COSR, $ SINR ) IF( NRU.GT.0 ) - $ CALL DROT( NRU, U( 1, M-1 ), 1, U( 1, M ), 1, COSL, SINL ) + $ CALL DROT( NRU, U( 1, M-1 ), 1, U( 1, M ), 1, COSL, + $ SINL ) IF( NCC.GT.0 ) $ CALL DROT( NCC, C( M-1, 1 ), LDC, C( M, 1 ), LDC, COSL, $ SINL ) @@ -629,7 +636,8 @@ SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, CALL DLARTG( D( I )*CS, E( I ), CS, SN, R ) IF( I.GT.LL ) $ E( I-1 ) = OLDSN*R - CALL DLARTG( OLDCS*R, D( I+1 )*SN, OLDCS, OLDSN, D( I ) ) + CALL DLARTG( OLDCS*R, D( I+1 )*SN, OLDCS, OLDSN, + $ D( I ) ) WORK( I-LL+1 ) = CS WORK( I-LL+1+NM1 ) = SN WORK( I-LL+1+NM12 ) = OLDCS @@ -645,10 +653,12 @@ SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, $ CALL DLASR( 'L', 'V', 'F', M-LL+1, NCVT, WORK( 1 ), $ WORK( N ), VT( LL, 1 ), LDVT ) IF( NRU.GT.0 ) - $ CALL DLASR( 'R', 'V', 'F', NRU, M-LL+1, WORK( NM12+1 ), + $ CALL DLASR( 'R', 'V', 'F', NRU, M-LL+1, + $ WORK( NM12+1 ), $ WORK( NM13+1 ), U( 1, LL ), LDU ) IF( NCC.GT.0 ) - $ CALL DLASR( 'L', 'V', 'F', M-LL+1, NCC, WORK( NM12+1 ), + $ CALL DLASR( 'L', 'V', 'F', M-LL+1, NCC, + $ WORK( NM12+1 ), $ WORK( NM13+1 ), C( LL, 1 ), LDC ) * * Test convergence @@ -667,7 +677,8 @@ SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, CALL DLARTG( D( I )*CS, E( I-1 ), CS, SN, R ) IF( I.LT.M ) $ E( I ) = OLDSN*R - CALL DLARTG( OLDCS*R, D( I-1 )*SN, OLDCS, OLDSN, D( I ) ) + CALL DLARTG( OLDCS*R, D( I-1 )*SN, OLDCS, OLDSN, + $ D( I ) ) WORK( I-LL ) = CS WORK( I-LL+NM1 ) = -SN WORK( I-LL+NM12 ) = OLDCS @@ -680,7 +691,8 @@ SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, * Update singular vectors * IF( NCVT.GT.0 ) - $ CALL DLASR( 'L', 'V', 'B', M-LL+1, NCVT, WORK( NM12+1 ), + $ CALL DLASR( 'L', 'V', 'B', M-LL+1, NCVT, + $ WORK( NM12+1 ), $ WORK( NM13+1 ), VT( LL, 1 ), LDVT ) IF( NRU.GT.0 ) $ CALL DLASR( 'R', 'V', 'B', NRU, M-LL+1, WORK( 1 ), @@ -735,10 +747,12 @@ SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, $ CALL DLASR( 'L', 'V', 'F', M-LL+1, NCVT, WORK( 1 ), $ WORK( N ), VT( LL, 1 ), LDVT ) IF( NRU.GT.0 ) - $ CALL DLASR( 'R', 'V', 'F', NRU, M-LL+1, WORK( NM12+1 ), + $ CALL DLASR( 'R', 'V', 'F', NRU, M-LL+1, + $ WORK( NM12+1 ), $ WORK( NM13+1 ), U( 1, LL ), LDU ) IF( NCC.GT.0 ) - $ CALL DLASR( 'L', 'V', 'F', M-LL+1, NCC, WORK( NM12+1 ), + $ CALL DLASR( 'L', 'V', 'F', M-LL+1, NCC, + $ WORK( NM12+1 ), $ WORK( NM13+1 ), C( LL, 1 ), LDC ) * * Test convergence @@ -785,7 +799,8 @@ SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, * Update singular vectors if desired * IF( NCVT.GT.0 ) - $ CALL DLASR( 'L', 'V', 'B', M-LL+1, NCVT, WORK( NM12+1 ), + $ CALL DLASR( 'L', 'V', 'B', M-LL+1, NCVT, + $ WORK( NM12+1 ), $ WORK( NM13+1 ), VT( LL, 1 ), LDVT ) IF( NRU.GT.0 ) $ CALL DLASR( 'R', 'V', 'B', NRU, M-LL+1, WORK( 1 ), @@ -804,6 +819,12 @@ SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, * 160 CONTINUE DO 170 I = 1, N + IF( D( I ).EQ.ZERO ) THEN +* +* Avoid -ZERO +* + D( I ) = ZERO + END IF IF( D( I ).LT.ZERO ) THEN D( I ) = -D( I ) * @@ -841,7 +862,8 @@ SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, IF( NRU.GT.0 ) $ CALL DSWAP( NRU, U( 1, ISUB ), 1, U( 1, N+1-I ), 1 ) IF( NCC.GT.0 ) - $ CALL DSWAP( NCC, C( ISUB, 1 ), LDC, C( N+1-I, 1 ), LDC ) + $ CALL DSWAP( NCC, C( ISUB, 1 ), LDC, C( N+1-I, 1 ), + $ LDC ) END IF 190 CONTINUE GO TO 220 diff --git a/lapack-netlib/SRC/sbdsqr.f b/lapack-netlib/SRC/sbdsqr.f index 880f0607b7..866b132950 100644 --- a/lapack-netlib/SRC/sbdsqr.f +++ b/lapack-netlib/SRC/sbdsqr.f @@ -5,7 +5,6 @@ * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> \htmlonly *> Download SBDSQR + dependencies *> *> [TGZ] @@ -13,7 +12,6 @@ *> [ZIP] *> *> [TXT] -*> \endhtmlonly * * Definition: * =========== @@ -166,7 +164,9 @@ *> *> \param[out] WORK *> \verbatim -*> WORK is REAL array, dimension (4*N) +*> WORK is REAL array, dimension (LWORK) +*> LWORK = 4*N, if NCVT = NRU = NCC = 0, and +*> LWORK = 4*(N-1), otherwise *> \endverbatim *> *> \param[out] INFO @@ -181,7 +181,7 @@ *> iterations (in inner while loop) *> = 3, termination criterion of outer while loop not met *> (program created more than N unreduced blocks) -*> else NCVT = NRU = NCC = 0, +*> else, *> the algorithm did not converge; D and E contain the *> elements of a bidiagonal matrix which is orthogonally *> similar to the input matrix B; if INFO = i, i @@ -232,11 +232,12 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup auxOTHERcomputational +*> \ingroup bdsqr * * ===================================================================== SUBROUTINE SBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, $ LDU, C, LDC, WORK, INFO ) + IMPLICIT NONE * * -- LAPACK computational routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -286,7 +287,8 @@ SUBROUTINE SBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, EXTERNAL LSAME, SLAMCH * .. * .. External Subroutines .. - EXTERNAL SLARTG, SLAS2, SLASQ1, SLASR, SLASV2, SROT, + EXTERNAL SLARTG, SLAS2, SLASQ1, SLASR, SLASV2, + $ SROT, $ SSCAL, SSWAP, XERBLA * .. * .. Intrinsic Functions .. @@ -335,9 +337,10 @@ SUBROUTINE SBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, IF( .NOT.ROTATE ) THEN CALL SLASQ1( N, D, E, WORK, INFO ) * -* If INFO equals 2, dqds didn't finish, try to finish +* If the dqds algorithm failed, try to finish with +* the standard QR algorithm * - IF( INFO .NE. 2 ) RETURN + IF( INFO .EQ. 0 ) RETURN INFO = 0 END IF * @@ -367,10 +370,12 @@ SUBROUTINE SBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, * Update singular vectors if desired * IF( NRU.GT.0 ) - $ CALL SLASR( 'R', 'V', 'F', NRU, N, WORK( 1 ), WORK( N ), U, + $ CALL SLASR( 'R', 'V', 'F', NRU, N, WORK( 1 ), WORK( N ), + $ U, $ LDU ) IF( NCC.GT.0 ) - $ CALL SLASR( 'L', 'V', 'F', N, NCC, WORK( 1 ), WORK( N ), C, + $ CALL SLASR( 'L', 'V', 'F', N, NCC, WORK( 1 ), WORK( N ), + $ C, $ LDC ) END IF * @@ -407,12 +412,13 @@ SUBROUTINE SBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, 40 CONTINUE 50 CONTINUE SMINOA = SMINOA / SQRT( REAL( N ) ) - THRESH = MAX( TOL*SMINOA, MAXITR*(N*(N*UNFL)) ) + THRESH = MAX( TOL*SMINOA, MAXITR*(REAL( N )*(REAL( N )*UNFL)) ) ELSE * * Absolute accuracy desired * - THRESH = MAX( ABS( TOL )*SMAX, MAXITR*(N*(N*UNFL)) ) + THRESH = MAX( ABS( TOL )*SMAX, MAXITR*(REAL( N )* + $ (REAL( N )*UNFL)) ) END IF * * Prepare for main iteration loop for the singular values @@ -492,10 +498,12 @@ SUBROUTINE SBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, * Compute singular vectors, if desired * IF( NCVT.GT.0 ) - $ CALL SROT( NCVT, VT( M-1, 1 ), LDVT, VT( M, 1 ), LDVT, COSR, + $ CALL SROT( NCVT, VT( M-1, 1 ), LDVT, VT( M, 1 ), LDVT, + $ COSR, $ SINR ) IF( NRU.GT.0 ) - $ CALL SROT( NRU, U( 1, M-1 ), 1, U( 1, M ), 1, COSL, SINL ) + $ CALL SROT( NRU, U( 1, M-1 ), 1, U( 1, M ), 1, COSL, + $ SINL ) IF( NCC.GT.0 ) $ CALL SROT( NCC, C( M-1, 1 ), LDC, C( M, 1 ), LDC, COSL, $ SINL ) @@ -584,7 +592,7 @@ SUBROUTINE SBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, * Compute shift. First, test if shifting would ruin relative * accuracy, and if so set the shift to zero. * - IF( TOL.GE.ZERO .AND. N*TOL*( SMIN / SMAX ).LE. + IF( TOL.GE.ZERO .AND. REAL( N )*TOL*( SMIN / SMAX ).LE. $ MAX( EPS, HNDRTH*TOL ) ) THEN * * Use a zero shift to avoid loss of relative accuracy @@ -628,7 +636,8 @@ SUBROUTINE SBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, CALL SLARTG( D( I )*CS, E( I ), CS, SN, R ) IF( I.GT.LL ) $ E( I-1 ) = OLDSN*R - CALL SLARTG( OLDCS*R, D( I+1 )*SN, OLDCS, OLDSN, D( I ) ) + CALL SLARTG( OLDCS*R, D( I+1 )*SN, OLDCS, OLDSN, + $ D( I ) ) WORK( I-LL+1 ) = CS WORK( I-LL+1+NM1 ) = SN WORK( I-LL+1+NM12 ) = OLDCS @@ -644,10 +653,12 @@ SUBROUTINE SBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, $ CALL SLASR( 'L', 'V', 'F', M-LL+1, NCVT, WORK( 1 ), $ WORK( N ), VT( LL, 1 ), LDVT ) IF( NRU.GT.0 ) - $ CALL SLASR( 'R', 'V', 'F', NRU, M-LL+1, WORK( NM12+1 ), + $ CALL SLASR( 'R', 'V', 'F', NRU, M-LL+1, + $ WORK( NM12+1 ), $ WORK( NM13+1 ), U( 1, LL ), LDU ) IF( NCC.GT.0 ) - $ CALL SLASR( 'L', 'V', 'F', M-LL+1, NCC, WORK( NM12+1 ), + $ CALL SLASR( 'L', 'V', 'F', M-LL+1, NCC, + $ WORK( NM12+1 ), $ WORK( NM13+1 ), C( LL, 1 ), LDC ) * * Test convergence @@ -666,7 +677,8 @@ SUBROUTINE SBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, CALL SLARTG( D( I )*CS, E( I-1 ), CS, SN, R ) IF( I.LT.M ) $ E( I ) = OLDSN*R - CALL SLARTG( OLDCS*R, D( I-1 )*SN, OLDCS, OLDSN, D( I ) ) + CALL SLARTG( OLDCS*R, D( I-1 )*SN, OLDCS, OLDSN, + $ D( I ) ) WORK( I-LL ) = CS WORK( I-LL+NM1 ) = -SN WORK( I-LL+NM12 ) = OLDCS @@ -679,7 +691,8 @@ SUBROUTINE SBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, * Update singular vectors * IF( NCVT.GT.0 ) - $ CALL SLASR( 'L', 'V', 'B', M-LL+1, NCVT, WORK( NM12+1 ), + $ CALL SLASR( 'L', 'V', 'B', M-LL+1, NCVT, + $ WORK( NM12+1 ), $ WORK( NM13+1 ), VT( LL, 1 ), LDVT ) IF( NRU.GT.0 ) $ CALL SLASR( 'R', 'V', 'B', NRU, M-LL+1, WORK( 1 ), @@ -734,10 +747,12 @@ SUBROUTINE SBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, $ CALL SLASR( 'L', 'V', 'F', M-LL+1, NCVT, WORK( 1 ), $ WORK( N ), VT( LL, 1 ), LDVT ) IF( NRU.GT.0 ) - $ CALL SLASR( 'R', 'V', 'F', NRU, M-LL+1, WORK( NM12+1 ), + $ CALL SLASR( 'R', 'V', 'F', NRU, M-LL+1, + $ WORK( NM12+1 ), $ WORK( NM13+1 ), U( 1, LL ), LDU ) IF( NCC.GT.0 ) - $ CALL SLASR( 'L', 'V', 'F', M-LL+1, NCC, WORK( NM12+1 ), + $ CALL SLASR( 'L', 'V', 'F', M-LL+1, NCC, + $ WORK( NM12+1 ), $ WORK( NM13+1 ), C( LL, 1 ), LDC ) * * Test convergence @@ -784,7 +799,8 @@ SUBROUTINE SBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, * Update singular vectors if desired * IF( NCVT.GT.0 ) - $ CALL SLASR( 'L', 'V', 'B', M-LL+1, NCVT, WORK( NM12+1 ), + $ CALL SLASR( 'L', 'V', 'B', M-LL+1, NCVT, + $ WORK( NM12+1 ), $ WORK( NM13+1 ), VT( LL, 1 ), LDVT ) IF( NRU.GT.0 ) $ CALL SLASR( 'R', 'V', 'B', NRU, M-LL+1, WORK( 1 ), @@ -803,7 +819,14 @@ SUBROUTINE SBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, * 160 CONTINUE DO 170 I = 1, N + IF( D( I ).EQ.ZERO ) THEN +* +* Avoid -ZERO +* + D( I ) = ZERO + END IF IF( D( I ).LT.ZERO ) THEN + D( I ) = -D( I ) * * Change sign of singular vectors, if desired @@ -840,7 +863,8 @@ SUBROUTINE SBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, IF( NRU.GT.0 ) $ CALL SSWAP( NRU, U( 1, ISUB ), 1, U( 1, N+1-I ), 1 ) IF( NCC.GT.0 ) - $ CALL SSWAP( NCC, C( ISUB, 1 ), LDC, C( N+1-I, 1 ), LDC ) + $ CALL SSWAP( NCC, C( ISUB, 1 ), LDC, C( N+1-I, 1 ), + $ LDC ) END IF 190 CONTINUE GO TO 220 diff --git a/lapack-netlib/SRC/zbdsqr.f b/lapack-netlib/SRC/zbdsqr.f index 865bb9dd59..ec67efeb02 100644 --- a/lapack-netlib/SRC/zbdsqr.f +++ b/lapack-netlib/SRC/zbdsqr.f @@ -5,7 +5,6 @@ * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> \htmlonly *> Download ZBDSQR + dependencies *> *> [TGZ] @@ -13,7 +12,6 @@ *> [ZIP] *> *> [TXT] -*> \endhtmlonly * * Definition: * =========== @@ -166,7 +164,9 @@ *> *> \param[out] RWORK *> \verbatim -*> RWORK is DOUBLE PRECISION array, dimension (4*N) +*> RWORK is DOUBLE PRECISION array, dimension (LRWORK) +*> LRWORK = 4*N, if NCVT = NRU = NCC = 0, and +*> LRWORK = 4*(N-1), otherwise *> \endverbatim *> *> \param[out] INFO @@ -230,6 +230,7 @@ * ===================================================================== SUBROUTINE ZBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, $ LDU, C, LDC, RWORK, INFO ) + IMPLICIT NONE * * -- LAPACK computational routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -279,7 +280,8 @@ SUBROUTINE ZBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, EXTERNAL LSAME, DLAMCH * .. * .. External Subroutines .. - EXTERNAL DLARTG, DLAS2, DLASQ1, DLASV2, XERBLA, ZDROT, + EXTERNAL DLARTG, DLAS2, DLASQ1, DLASV2, XERBLA, + $ ZDROT, $ ZDSCAL, ZLASR, ZSWAP * .. * .. Intrinsic Functions .. @@ -328,9 +330,10 @@ SUBROUTINE ZBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, IF( .NOT.ROTATE ) THEN CALL DLASQ1( N, D, E, RWORK, INFO ) * -* If INFO equals 2, dqds didn't finish, try to finish +* If the dqds algorithm failed, try to finish with +* the standard QR algorithm * - IF( INFO .NE. 2 ) RETURN + IF( INFO .EQ. 0 ) RETURN INFO = 0 END IF * @@ -360,10 +363,12 @@ SUBROUTINE ZBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, * Update singular vectors if desired * IF( NRU.GT.0 ) - $ CALL ZLASR( 'R', 'V', 'F', NRU, N, RWORK( 1 ), RWORK( N ), + $ CALL ZLASR( 'R', 'V', 'F', NRU, N, RWORK( 1 ), + $ RWORK( N ), $ U, LDU ) IF( NCC.GT.0 ) - $ CALL ZLASR( 'L', 'V', 'F', N, NCC, RWORK( 1 ), RWORK( N ), + $ CALL ZLASR( 'L', 'V', 'F', N, NCC, RWORK( 1 ), + $ RWORK( N ), $ C, LDC ) END IF * @@ -487,7 +492,8 @@ SUBROUTINE ZBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, $ CALL ZDROT( NCVT, VT( M-1, 1 ), LDVT, VT( M, 1 ), LDVT, $ COSR, SINR ) IF( NRU.GT.0 ) - $ CALL ZDROT( NRU, U( 1, M-1 ), 1, U( 1, M ), 1, COSL, SINL ) + $ CALL ZDROT( NRU, U( 1, M-1 ), 1, U( 1, M ), 1, COSL, + $ SINL ) IF( NCC.GT.0 ) $ CALL ZDROT( NCC, C( M-1, 1 ), LDC, C( M, 1 ), LDC, COSL, $ SINL ) @@ -620,7 +626,8 @@ SUBROUTINE ZBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, CALL DLARTG( D( I )*CS, E( I ), CS, SN, R ) IF( I.GT.LL ) $ E( I-1 ) = OLDSN*R - CALL DLARTG( OLDCS*R, D( I+1 )*SN, OLDCS, OLDSN, D( I ) ) + CALL DLARTG( OLDCS*R, D( I+1 )*SN, OLDCS, OLDSN, + $ D( I ) ) RWORK( I-LL+1 ) = CS RWORK( I-LL+1+NM1 ) = SN RWORK( I-LL+1+NM12 ) = OLDCS @@ -636,10 +643,12 @@ SUBROUTINE ZBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, $ CALL ZLASR( 'L', 'V', 'F', M-LL+1, NCVT, RWORK( 1 ), $ RWORK( N ), VT( LL, 1 ), LDVT ) IF( NRU.GT.0 ) - $ CALL ZLASR( 'R', 'V', 'F', NRU, M-LL+1, RWORK( NM12+1 ), + $ CALL ZLASR( 'R', 'V', 'F', NRU, M-LL+1, + $ RWORK( NM12+1 ), $ RWORK( NM13+1 ), U( 1, LL ), LDU ) IF( NCC.GT.0 ) - $ CALL ZLASR( 'L', 'V', 'F', M-LL+1, NCC, RWORK( NM12+1 ), + $ CALL ZLASR( 'L', 'V', 'F', M-LL+1, NCC, + $ RWORK( NM12+1 ), $ RWORK( NM13+1 ), C( LL, 1 ), LDC ) * * Test convergence @@ -658,7 +667,8 @@ SUBROUTINE ZBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, CALL DLARTG( D( I )*CS, E( I-1 ), CS, SN, R ) IF( I.LT.M ) $ E( I ) = OLDSN*R - CALL DLARTG( OLDCS*R, D( I-1 )*SN, OLDCS, OLDSN, D( I ) ) + CALL DLARTG( OLDCS*R, D( I-1 )*SN, OLDCS, OLDSN, + $ D( I ) ) RWORK( I-LL ) = CS RWORK( I-LL+NM1 ) = -SN RWORK( I-LL+NM12 ) = OLDCS @@ -671,7 +681,8 @@ SUBROUTINE ZBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, * Update singular vectors * IF( NCVT.GT.0 ) - $ CALL ZLASR( 'L', 'V', 'B', M-LL+1, NCVT, RWORK( NM12+1 ), + $ CALL ZLASR( 'L', 'V', 'B', M-LL+1, NCVT, + $ RWORK( NM12+1 ), $ RWORK( NM13+1 ), VT( LL, 1 ), LDVT ) IF( NRU.GT.0 ) $ CALL ZLASR( 'R', 'V', 'B', NRU, M-LL+1, RWORK( 1 ), @@ -726,10 +737,12 @@ SUBROUTINE ZBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, $ CALL ZLASR( 'L', 'V', 'F', M-LL+1, NCVT, RWORK( 1 ), $ RWORK( N ), VT( LL, 1 ), LDVT ) IF( NRU.GT.0 ) - $ CALL ZLASR( 'R', 'V', 'F', NRU, M-LL+1, RWORK( NM12+1 ), + $ CALL ZLASR( 'R', 'V', 'F', NRU, M-LL+1, + $ RWORK( NM12+1 ), $ RWORK( NM13+1 ), U( 1, LL ), LDU ) IF( NCC.GT.0 ) - $ CALL ZLASR( 'L', 'V', 'F', M-LL+1, NCC, RWORK( NM12+1 ), + $ CALL ZLASR( 'L', 'V', 'F', M-LL+1, NCC, + $ RWORK( NM12+1 ), $ RWORK( NM13+1 ), C( LL, 1 ), LDC ) * * Test convergence @@ -776,7 +789,8 @@ SUBROUTINE ZBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, * Update singular vectors if desired * IF( NCVT.GT.0 ) - $ CALL ZLASR( 'L', 'V', 'B', M-LL+1, NCVT, RWORK( NM12+1 ), + $ CALL ZLASR( 'L', 'V', 'B', M-LL+1, NCVT, + $ RWORK( NM12+1 ), $ RWORK( NM13+1 ), VT( LL, 1 ), LDVT ) IF( NRU.GT.0 ) $ CALL ZLASR( 'R', 'V', 'B', NRU, M-LL+1, RWORK( 1 ), @@ -795,6 +809,12 @@ SUBROUTINE ZBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, * 160 CONTINUE DO 170 I = 1, N + IF( D( I ).EQ.ZERO ) THEN +* +* Avoid -ZERO +* + D( I ) = ZERO + END IF IF( D( I ).LT.ZERO ) THEN D( I ) = -D( I ) * @@ -832,7 +852,8 @@ SUBROUTINE ZBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, IF( NRU.GT.0 ) $ CALL ZSWAP( NRU, U( 1, ISUB ), 1, U( 1, N+1-I ), 1 ) IF( NCC.GT.0 ) - $ CALL ZSWAP( NCC, C( ISUB, 1 ), LDC, C( N+1-I, 1 ), LDC ) + $ CALL ZSWAP( NCC, C( ISUB, 1 ), LDC, C( N+1-I, 1 ), + $ LDC ) END IF 190 CONTINUE GO TO 220