summaryrefslogtreecommitdiff
path: root/src/fortran/lapack/zbdsqr.f
diff options
context:
space:
mode:
Diffstat (limited to 'src/fortran/lapack/zbdsqr.f')
-rw-r--r--src/fortran/lapack/zbdsqr.f742
1 files changed, 742 insertions, 0 deletions
diff --git a/src/fortran/lapack/zbdsqr.f b/src/fortran/lapack/zbdsqr.f
new file mode 100644
index 0000000..f9086be
--- /dev/null
+++ b/src/fortran/lapack/zbdsqr.f
@@ -0,0 +1,742 @@
+ SUBROUTINE ZBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
+ $ LDU, C, LDC, RWORK, INFO )
+*
+* -- LAPACK routine (version 3.1) --
+* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+* November 2006
+*
+* .. Scalar Arguments ..
+ CHARACTER UPLO
+ INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU
+* ..
+* .. Array Arguments ..
+ DOUBLE PRECISION D( * ), E( * ), RWORK( * )
+ COMPLEX*16 C( LDC, * ), U( LDU, * ), VT( LDVT, * )
+* ..
+*
+* Purpose
+* =======
+*
+* ZBDSQR computes the singular values and, optionally, the right and/or
+* left singular vectors from the singular value decomposition (SVD) of
+* a real N-by-N (upper or lower) bidiagonal matrix B using the implicit
+* zero-shift QR algorithm. The SVD of B has the form
+*
+* B = Q * S * P**H
+*
+* where S is the diagonal matrix of singular values, Q is an orthogonal
+* matrix of left singular vectors, and P is an orthogonal matrix of
+* right singular vectors. If left singular vectors are requested, this
+* subroutine actually returns U*Q instead of Q, and, if right singular
+* vectors are requested, this subroutine returns P**H*VT instead of
+* P**H, for given complex input matrices U and VT. When U and VT are
+* the unitary matrices that reduce a general matrix A to bidiagonal
+* form: A = U*B*VT, as computed by ZGEBRD, then
+*
+* A = (U*Q) * S * (P**H*VT)
+*
+* is the SVD of A. Optionally, the subroutine may also compute Q**H*C
+* for a given complex input matrix C.
+*
+* See "Computing Small Singular Values of Bidiagonal Matrices With
+* Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
+* LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,
+* no. 5, pp. 873-912, Sept 1990) and
+* "Accurate singular values and differential qd algorithms," by
+* B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics
+* Department, University of California at Berkeley, July 1992
+* for a detailed description of the algorithm.
+*
+* Arguments
+* =========
+*
+* UPLO (input) CHARACTER*1
+* = 'U': B is upper bidiagonal;
+* = 'L': B is lower bidiagonal.
+*
+* N (input) INTEGER
+* The order of the matrix B. N >= 0.
+*
+* NCVT (input) INTEGER
+* The number of columns of the matrix VT. NCVT >= 0.
+*
+* NRU (input) INTEGER
+* The number of rows of the matrix U. NRU >= 0.
+*
+* NCC (input) INTEGER
+* The number of columns of the matrix C. NCC >= 0.
+*
+* D (input/output) DOUBLE PRECISION array, dimension (N)
+* On entry, the n diagonal elements of the bidiagonal matrix B.
+* On exit, if INFO=0, the singular values of B in decreasing
+* order.
+*
+* E (input/output) DOUBLE PRECISION array, dimension (N-1)
+* On entry, the N-1 offdiagonal elements of the bidiagonal
+* matrix B.
+* On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E
+* will contain the diagonal and superdiagonal elements of a
+* bidiagonal matrix orthogonally equivalent to the one given
+* as input.
+*
+* VT (input/output) COMPLEX*16 array, dimension (LDVT, NCVT)
+* On entry, an N-by-NCVT matrix VT.
+* On exit, VT is overwritten by P**H * VT.
+* Not referenced if NCVT = 0.
+*
+* LDVT (input) INTEGER
+* The leading dimension of the array VT.
+* LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0.
+*
+* U (input/output) COMPLEX*16 array, dimension (LDU, N)
+* On entry, an NRU-by-N matrix U.
+* On exit, U is overwritten by U * Q.
+* Not referenced if NRU = 0.
+*
+* LDU (input) INTEGER
+* The leading dimension of the array U. LDU >= max(1,NRU).
+*
+* C (input/output) COMPLEX*16 array, dimension (LDC, NCC)
+* On entry, an N-by-NCC matrix C.
+* On exit, C is overwritten by Q**H * C.
+* Not referenced if NCC = 0.
+*
+* LDC (input) INTEGER
+* The leading dimension of the array C.
+* LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.
+*
+* RWORK (workspace) DOUBLE PRECISION array, dimension (2*N)
+* if NCVT = NRU = NCC = 0, (max(1, 4*N-4)) otherwise
+*
+* INFO (output) INTEGER
+* = 0: successful exit
+* < 0: If INFO = -i, the i-th argument had an illegal value
+* > 0: 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
+* elements of E have not converged to zero.
+*
+* Internal Parameters
+* ===================
+*
+* TOLMUL DOUBLE PRECISION, default = max(10,min(100,EPS**(-1/8)))
+* TOLMUL controls the convergence criterion of the QR loop.
+* If it is positive, TOLMUL*EPS is the desired relative
+* precision in the computed singular values.
+* If it is negative, abs(TOLMUL*EPS*sigma_max) is the
+* desired absolute accuracy in the computed singular
+* values (corresponds to relative accuracy
+* abs(TOLMUL*EPS) in the largest singular value.
+* abs(TOLMUL) should be between 1 and 1/EPS, and preferably
+* between 10 (for fast convergence) and .1/EPS
+* (for there to be some accuracy in the results).
+* Default is to lose at either one eighth or 2 of the
+* available decimal digits in each computed singular value
+* (whichever is smaller).
+*
+* MAXITR INTEGER, default = 6
+* MAXITR controls the maximum number of passes of the
+* algorithm through its inner loop. The algorithms stops
+* (and so fails to converge) if the number of passes
+* through the inner loop exceeds MAXITR*N**2.
+*
+* =====================================================================
+*
+* .. Parameters ..
+ DOUBLE PRECISION ZERO
+ PARAMETER ( ZERO = 0.0D0 )
+ DOUBLE PRECISION ONE
+ PARAMETER ( ONE = 1.0D0 )
+ DOUBLE PRECISION NEGONE
+ PARAMETER ( NEGONE = -1.0D0 )
+ DOUBLE PRECISION HNDRTH
+ PARAMETER ( HNDRTH = 0.01D0 )
+ DOUBLE PRECISION TEN
+ PARAMETER ( TEN = 10.0D0 )
+ DOUBLE PRECISION HNDRD
+ PARAMETER ( HNDRD = 100.0D0 )
+ DOUBLE PRECISION MEIGTH
+ PARAMETER ( MEIGTH = -0.125D0 )
+ INTEGER MAXITR
+ PARAMETER ( MAXITR = 6 )
+* ..
+* .. Local Scalars ..
+ LOGICAL LOWER, ROTATE
+ INTEGER I, IDIR, ISUB, ITER, J, LL, LLL, M, MAXIT, NM1,
+ $ NM12, NM13, OLDLL, OLDM
+ DOUBLE PRECISION ABSE, ABSS, COSL, COSR, CS, EPS, F, G, H, MU,
+ $ OLDCS, OLDSN, R, SHIFT, SIGMN, SIGMX, SINL,
+ $ SINR, SLL, SMAX, SMIN, SMINL, SMINOA,
+ $ SN, THRESH, TOL, TOLMUL, UNFL
+* ..
+* .. External Functions ..
+ LOGICAL LSAME
+ DOUBLE PRECISION DLAMCH
+ EXTERNAL LSAME, DLAMCH
+* ..
+* .. External Subroutines ..
+ EXTERNAL DLARTG, DLAS2, DLASQ1, DLASV2, XERBLA, ZDROT,
+ $ ZDSCAL, ZLASR, ZSWAP
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC ABS, DBLE, MAX, MIN, SIGN, SQRT
+* ..
+* .. Executable Statements ..
+*
+* Test the input parameters.
+*
+ INFO = 0
+ LOWER = LSAME( UPLO, 'L' )
+ IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LOWER ) THEN
+ INFO = -1
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -2
+ ELSE IF( NCVT.LT.0 ) THEN
+ INFO = -3
+ ELSE IF( NRU.LT.0 ) THEN
+ INFO = -4
+ ELSE IF( NCC.LT.0 ) THEN
+ INFO = -5
+ ELSE IF( ( NCVT.EQ.0 .AND. LDVT.LT.1 ) .OR.
+ $ ( NCVT.GT.0 .AND. LDVT.LT.MAX( 1, N ) ) ) THEN
+ INFO = -9
+ ELSE IF( LDU.LT.MAX( 1, NRU ) ) THEN
+ INFO = -11
+ ELSE IF( ( NCC.EQ.0 .AND. LDC.LT.1 ) .OR.
+ $ ( NCC.GT.0 .AND. LDC.LT.MAX( 1, N ) ) ) THEN
+ INFO = -13
+ END IF
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'ZBDSQR', -INFO )
+ RETURN
+ END IF
+ IF( N.EQ.0 )
+ $ RETURN
+ IF( N.EQ.1 )
+ $ GO TO 160
+*
+* ROTATE is true if any singular vectors desired, false otherwise
+*
+ ROTATE = ( NCVT.GT.0 ) .OR. ( NRU.GT.0 ) .OR. ( NCC.GT.0 )
+*
+* If no singular vectors desired, use qd algorithm
+*
+ IF( .NOT.ROTATE ) THEN
+ CALL DLASQ1( N, D, E, RWORK, INFO )
+ RETURN
+ END IF
+*
+ NM1 = N - 1
+ NM12 = NM1 + NM1
+ NM13 = NM12 + NM1
+ IDIR = 0
+*
+* Get machine constants
+*
+ EPS = DLAMCH( 'Epsilon' )
+ UNFL = DLAMCH( 'Safe minimum' )
+*
+* If matrix lower bidiagonal, rotate to be upper bidiagonal
+* by applying Givens rotations on the left
+*
+ IF( LOWER ) THEN
+ DO 10 I = 1, N - 1
+ CALL DLARTG( D( I ), E( I ), CS, SN, R )
+ D( I ) = R
+ E( I ) = SN*D( I+1 )
+ D( I+1 ) = CS*D( I+1 )
+ RWORK( I ) = CS
+ RWORK( NM1+I ) = SN
+ 10 CONTINUE
+*
+* Update singular vectors if desired
+*
+ IF( NRU.GT.0 )
+ $ 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 ),
+ $ C, LDC )
+ END IF
+*
+* Compute singular values to relative accuracy TOL
+* (By setting TOL to be negative, algorithm will compute
+* singular values to absolute accuracy ABS(TOL)*norm(input matrix))
+*
+ TOLMUL = MAX( TEN, MIN( HNDRD, EPS**MEIGTH ) )
+ TOL = TOLMUL*EPS
+*
+* Compute approximate maximum, minimum singular values
+*
+ SMAX = ZERO
+ DO 20 I = 1, N
+ SMAX = MAX( SMAX, ABS( D( I ) ) )
+ 20 CONTINUE
+ DO 30 I = 1, N - 1
+ SMAX = MAX( SMAX, ABS( E( I ) ) )
+ 30 CONTINUE
+ SMINL = ZERO
+ IF( TOL.GE.ZERO ) THEN
+*
+* Relative accuracy desired
+*
+ SMINOA = ABS( D( 1 ) )
+ IF( SMINOA.EQ.ZERO )
+ $ GO TO 50
+ MU = SMINOA
+ DO 40 I = 2, N
+ MU = ABS( D( I ) )*( MU / ( MU+ABS( E( I-1 ) ) ) )
+ SMINOA = MIN( SMINOA, MU )
+ IF( SMINOA.EQ.ZERO )
+ $ GO TO 50
+ 40 CONTINUE
+ 50 CONTINUE
+ SMINOA = SMINOA / SQRT( DBLE( N ) )
+ THRESH = MAX( TOL*SMINOA, MAXITR*N*N*UNFL )
+ ELSE
+*
+* Absolute accuracy desired
+*
+ THRESH = MAX( ABS( TOL )*SMAX, MAXITR*N*N*UNFL )
+ END IF
+*
+* Prepare for main iteration loop for the singular values
+* (MAXIT is the maximum number of passes through the inner
+* loop permitted before nonconvergence signalled.)
+*
+ MAXIT = MAXITR*N*N
+ ITER = 0
+ OLDLL = -1
+ OLDM = -1
+*
+* M points to last element of unconverged part of matrix
+*
+ M = N
+*
+* Begin main iteration loop
+*
+ 60 CONTINUE
+*
+* Check for convergence or exceeding iteration count
+*
+ IF( M.LE.1 )
+ $ GO TO 160
+ IF( ITER.GT.MAXIT )
+ $ GO TO 200
+*
+* Find diagonal block of matrix to work on
+*
+ IF( TOL.LT.ZERO .AND. ABS( D( M ) ).LE.THRESH )
+ $ D( M ) = ZERO
+ SMAX = ABS( D( M ) )
+ SMIN = SMAX
+ DO 70 LLL = 1, M - 1
+ LL = M - LLL
+ ABSS = ABS( D( LL ) )
+ ABSE = ABS( E( LL ) )
+ IF( TOL.LT.ZERO .AND. ABSS.LE.THRESH )
+ $ D( LL ) = ZERO
+ IF( ABSE.LE.THRESH )
+ $ GO TO 80
+ SMIN = MIN( SMIN, ABSS )
+ SMAX = MAX( SMAX, ABSS, ABSE )
+ 70 CONTINUE
+ LL = 0
+ GO TO 90
+ 80 CONTINUE
+ E( LL ) = ZERO
+*
+* Matrix splits since E(LL) = 0
+*
+ IF( LL.EQ.M-1 ) THEN
+*
+* Convergence of bottom singular value, return to top of loop
+*
+ M = M - 1
+ GO TO 60
+ END IF
+ 90 CONTINUE
+ LL = LL + 1
+*
+* E(LL) through E(M-1) are nonzero, E(LL-1) is zero
+*
+ IF( LL.EQ.M-1 ) THEN
+*
+* 2 by 2 block, handle separately
+*
+ CALL DLASV2( D( M-1 ), E( M-1 ), D( M ), SIGMN, SIGMX, SINR,
+ $ COSR, SINL, COSL )
+ D( M-1 ) = SIGMX
+ E( M-1 ) = ZERO
+ D( M ) = SIGMN
+*
+* Compute singular vectors, if desired
+*
+ IF( NCVT.GT.0 )
+ $ 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 )
+ IF( NCC.GT.0 )
+ $ CALL ZDROT( NCC, C( M-1, 1 ), LDC, C( M, 1 ), LDC, COSL,
+ $ SINL )
+ M = M - 2
+ GO TO 60
+ END IF
+*
+* If working on new submatrix, choose shift direction
+* (from larger end diagonal element towards smaller)
+*
+ IF( LL.GT.OLDM .OR. M.LT.OLDLL ) THEN
+ IF( ABS( D( LL ) ).GE.ABS( D( M ) ) ) THEN
+*
+* Chase bulge from top (big end) to bottom (small end)
+*
+ IDIR = 1
+ ELSE
+*
+* Chase bulge from bottom (big end) to top (small end)
+*
+ IDIR = 2
+ END IF
+ END IF
+*
+* Apply convergence tests
+*
+ IF( IDIR.EQ.1 ) THEN
+*
+* Run convergence test in forward direction
+* First apply standard test to bottom of matrix
+*
+ IF( ABS( E( M-1 ) ).LE.ABS( TOL )*ABS( D( M ) ) .OR.
+ $ ( TOL.LT.ZERO .AND. ABS( E( M-1 ) ).LE.THRESH ) ) THEN
+ E( M-1 ) = ZERO
+ GO TO 60
+ END IF
+*
+ IF( TOL.GE.ZERO ) THEN
+*
+* If relative accuracy desired,
+* apply convergence criterion forward
+*
+ MU = ABS( D( LL ) )
+ SMINL = MU
+ DO 100 LLL = LL, M - 1
+ IF( ABS( E( LLL ) ).LE.TOL*MU ) THEN
+ E( LLL ) = ZERO
+ GO TO 60
+ END IF
+ MU = ABS( D( LLL+1 ) )*( MU / ( MU+ABS( E( LLL ) ) ) )
+ SMINL = MIN( SMINL, MU )
+ 100 CONTINUE
+ END IF
+*
+ ELSE
+*
+* Run convergence test in backward direction
+* First apply standard test to top of matrix
+*
+ IF( ABS( E( LL ) ).LE.ABS( TOL )*ABS( D( LL ) ) .OR.
+ $ ( TOL.LT.ZERO .AND. ABS( E( LL ) ).LE.THRESH ) ) THEN
+ E( LL ) = ZERO
+ GO TO 60
+ END IF
+*
+ IF( TOL.GE.ZERO ) THEN
+*
+* If relative accuracy desired,
+* apply convergence criterion backward
+*
+ MU = ABS( D( M ) )
+ SMINL = MU
+ DO 110 LLL = M - 1, LL, -1
+ IF( ABS( E( LLL ) ).LE.TOL*MU ) THEN
+ E( LLL ) = ZERO
+ GO TO 60
+ END IF
+ MU = ABS( D( LLL ) )*( MU / ( MU+ABS( E( LLL ) ) ) )
+ SMINL = MIN( SMINL, MU )
+ 110 CONTINUE
+ END IF
+ END IF
+ OLDLL = LL
+ OLDM = M
+*
+* 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*( SMINL / SMAX ).LE.
+ $ MAX( EPS, HNDRTH*TOL ) ) THEN
+*
+* Use a zero shift to avoid loss of relative accuracy
+*
+ SHIFT = ZERO
+ ELSE
+*
+* Compute the shift from 2-by-2 block at end of matrix
+*
+ IF( IDIR.EQ.1 ) THEN
+ SLL = ABS( D( LL ) )
+ CALL DLAS2( D( M-1 ), E( M-1 ), D( M ), SHIFT, R )
+ ELSE
+ SLL = ABS( D( M ) )
+ CALL DLAS2( D( LL ), E( LL ), D( LL+1 ), SHIFT, R )
+ END IF
+*
+* Test if shift negligible, and if so set to zero
+*
+ IF( SLL.GT.ZERO ) THEN
+ IF( ( SHIFT / SLL )**2.LT.EPS )
+ $ SHIFT = ZERO
+ END IF
+ END IF
+*
+* Increment iteration count
+*
+ ITER = ITER + M - LL
+*
+* If SHIFT = 0, do simplified QR iteration
+*
+ IF( SHIFT.EQ.ZERO ) THEN
+ IF( IDIR.EQ.1 ) THEN
+*
+* Chase bulge from top to bottom
+* Save cosines and sines for later singular vector updates
+*
+ CS = ONE
+ OLDCS = ONE
+ DO 120 I = LL, M - 1
+ 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 ) )
+ RWORK( I-LL+1 ) = CS
+ RWORK( I-LL+1+NM1 ) = SN
+ RWORK( I-LL+1+NM12 ) = OLDCS
+ RWORK( I-LL+1+NM13 ) = OLDSN
+ 120 CONTINUE
+ H = D( M )*CS
+ D( M ) = H*OLDCS
+ E( M-1 ) = H*OLDSN
+*
+* Update singular vectors
+*
+ IF( NCVT.GT.0 )
+ $ 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 ),
+ $ RWORK( NM13+1 ), U( 1, LL ), LDU )
+ IF( NCC.GT.0 )
+ $ CALL ZLASR( 'L', 'V', 'F', M-LL+1, NCC, RWORK( NM12+1 ),
+ $ RWORK( NM13+1 ), C( LL, 1 ), LDC )
+*
+* Test convergence
+*
+ IF( ABS( E( M-1 ) ).LE.THRESH )
+ $ E( M-1 ) = ZERO
+*
+ ELSE
+*
+* Chase bulge from bottom to top
+* Save cosines and sines for later singular vector updates
+*
+ CS = ONE
+ OLDCS = ONE
+ DO 130 I = M, LL + 1, -1
+ 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 ) )
+ RWORK( I-LL ) = CS
+ RWORK( I-LL+NM1 ) = -SN
+ RWORK( I-LL+NM12 ) = OLDCS
+ RWORK( I-LL+NM13 ) = -OLDSN
+ 130 CONTINUE
+ H = D( LL )*CS
+ D( LL ) = H*OLDCS
+ E( LL ) = H*OLDSN
+*
+* Update singular vectors
+*
+ IF( NCVT.GT.0 )
+ $ 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 ),
+ $ RWORK( N ), U( 1, LL ), LDU )
+ IF( NCC.GT.0 )
+ $ CALL ZLASR( 'L', 'V', 'B', M-LL+1, NCC, RWORK( 1 ),
+ $ RWORK( N ), C( LL, 1 ), LDC )
+*
+* Test convergence
+*
+ IF( ABS( E( LL ) ).LE.THRESH )
+ $ E( LL ) = ZERO
+ END IF
+ ELSE
+*
+* Use nonzero shift
+*
+ IF( IDIR.EQ.1 ) THEN
+*
+* Chase bulge from top to bottom
+* Save cosines and sines for later singular vector updates
+*
+ F = ( ABS( D( LL ) )-SHIFT )*
+ $ ( SIGN( ONE, D( LL ) )+SHIFT / D( LL ) )
+ G = E( LL )
+ DO 140 I = LL, M - 1
+ CALL DLARTG( F, G, COSR, SINR, R )
+ IF( I.GT.LL )
+ $ E( I-1 ) = R
+ F = COSR*D( I ) + SINR*E( I )
+ E( I ) = COSR*E( I ) - SINR*D( I )
+ G = SINR*D( I+1 )
+ D( I+1 ) = COSR*D( I+1 )
+ CALL DLARTG( F, G, COSL, SINL, R )
+ D( I ) = R
+ F = COSL*E( I ) + SINL*D( I+1 )
+ D( I+1 ) = COSL*D( I+1 ) - SINL*E( I )
+ IF( I.LT.M-1 ) THEN
+ G = SINL*E( I+1 )
+ E( I+1 ) = COSL*E( I+1 )
+ END IF
+ RWORK( I-LL+1 ) = COSR
+ RWORK( I-LL+1+NM1 ) = SINR
+ RWORK( I-LL+1+NM12 ) = COSL
+ RWORK( I-LL+1+NM13 ) = SINL
+ 140 CONTINUE
+ E( M-1 ) = F
+*
+* Update singular vectors
+*
+ IF( NCVT.GT.0 )
+ $ 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 ),
+ $ RWORK( NM13+1 ), U( 1, LL ), LDU )
+ IF( NCC.GT.0 )
+ $ CALL ZLASR( 'L', 'V', 'F', M-LL+1, NCC, RWORK( NM12+1 ),
+ $ RWORK( NM13+1 ), C( LL, 1 ), LDC )
+*
+* Test convergence
+*
+ IF( ABS( E( M-1 ) ).LE.THRESH )
+ $ E( M-1 ) = ZERO
+*
+ ELSE
+*
+* Chase bulge from bottom to top
+* Save cosines and sines for later singular vector updates
+*
+ F = ( ABS( D( M ) )-SHIFT )*( SIGN( ONE, D( M ) )+SHIFT /
+ $ D( M ) )
+ G = E( M-1 )
+ DO 150 I = M, LL + 1, -1
+ CALL DLARTG( F, G, COSR, SINR, R )
+ IF( I.LT.M )
+ $ E( I ) = R
+ F = COSR*D( I ) + SINR*E( I-1 )
+ E( I-1 ) = COSR*E( I-1 ) - SINR*D( I )
+ G = SINR*D( I-1 )
+ D( I-1 ) = COSR*D( I-1 )
+ CALL DLARTG( F, G, COSL, SINL, R )
+ D( I ) = R
+ F = COSL*E( I-1 ) + SINL*D( I-1 )
+ D( I-1 ) = COSL*D( I-1 ) - SINL*E( I-1 )
+ IF( I.GT.LL+1 ) THEN
+ G = SINL*E( I-2 )
+ E( I-2 ) = COSL*E( I-2 )
+ END IF
+ RWORK( I-LL ) = COSR
+ RWORK( I-LL+NM1 ) = -SINR
+ RWORK( I-LL+NM12 ) = COSL
+ RWORK( I-LL+NM13 ) = -SINL
+ 150 CONTINUE
+ E( LL ) = F
+*
+* Test convergence
+*
+ IF( ABS( E( LL ) ).LE.THRESH )
+ $ E( LL ) = ZERO
+*
+* Update singular vectors if desired
+*
+ IF( NCVT.GT.0 )
+ $ 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 ),
+ $ RWORK( N ), U( 1, LL ), LDU )
+ IF( NCC.GT.0 )
+ $ CALL ZLASR( 'L', 'V', 'B', M-LL+1, NCC, RWORK( 1 ),
+ $ RWORK( N ), C( LL, 1 ), LDC )
+ END IF
+ END IF
+*
+* QR iteration finished, go back and check convergence
+*
+ GO TO 60
+*
+* All singular values converged, so make them positive
+*
+ 160 CONTINUE
+ DO 170 I = 1, N
+ IF( D( I ).LT.ZERO ) THEN
+ D( I ) = -D( I )
+*
+* Change sign of singular vectors, if desired
+*
+ IF( NCVT.GT.0 )
+ $ CALL ZDSCAL( NCVT, NEGONE, VT( I, 1 ), LDVT )
+ END IF
+ 170 CONTINUE
+*
+* Sort the singular values into decreasing order (insertion sort on
+* singular values, but only one transposition per singular vector)
+*
+ DO 190 I = 1, N - 1
+*
+* Scan for smallest D(I)
+*
+ ISUB = 1
+ SMIN = D( 1 )
+ DO 180 J = 2, N + 1 - I
+ IF( D( J ).LE.SMIN ) THEN
+ ISUB = J
+ SMIN = D( J )
+ END IF
+ 180 CONTINUE
+ IF( ISUB.NE.N+1-I ) THEN
+*
+* Swap singular values and vectors
+*
+ D( ISUB ) = D( N+1-I )
+ D( N+1-I ) = SMIN
+ IF( NCVT.GT.0 )
+ $ CALL ZSWAP( NCVT, VT( ISUB, 1 ), LDVT, VT( N+1-I, 1 ),
+ $ LDVT )
+ 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 )
+ END IF
+ 190 CONTINUE
+ GO TO 220
+*
+* Maximum number of iterations exceeded, failure to converge
+*
+ 200 CONTINUE
+ INFO = 0
+ DO 210 I = 1, N - 1
+ IF( E( I ).NE.ZERO )
+ $ INFO = INFO + 1
+ 210 CONTINUE
+ 220 CONTINUE
+ RETURN
+*
+* End of ZBDSQR
+*
+ END
generated by cgit (git 2.34.1) at 2024-12-20 18:20:26 +0000