From 8c8d2f518968ce7057eec6aa5cd5aec8faab861a Mon Sep 17 00:00:00 2001 From: jofret Date: Tue, 28 Apr 2009 07:17:00 +0000 Subject: Moving lapack to right place --- src/lib/lapack/zbdsqr.f | 742 ------------------------------------------------ 1 file changed, 742 deletions(-) delete mode 100644 src/lib/lapack/zbdsqr.f (limited to 'src/lib/lapack/zbdsqr.f') diff --git a/src/lib/lapack/zbdsqr.f b/src/lib/lapack/zbdsqr.f deleted file mode 100644 index f9086be5..00000000 --- a/src/lib/lapack/zbdsqr.f +++ /dev/null @@ -1,742 +0,0 @@ - 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 -- cgit