SUBROUTINE DLASDQ( UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT, \$ U, LDU, C, LDC, WORK, INFO ) * * -- LAPACK auxiliary 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, SQRE * .. * .. Array Arguments .. DOUBLE PRECISION C( LDC, * ), D( * ), E( * ), U( LDU, * ), \$ VT( LDVT, * ), WORK( * ) * .. * * Purpose * ======= * * DLASDQ computes the singular value decomposition (SVD) of a real * (upper or lower) bidiagonal matrix with diagonal D and offdiagonal * E, accumulating the transformations if desired. Letting B denote * the input bidiagonal matrix, the algorithm computes orthogonal * matrices Q and P such that B = Q * S * P' (P' denotes the transpose * of P). The singular values S are overwritten on D. * * The input matrix U is changed to U * Q if desired. * The input matrix VT is changed to P' * VT if desired. * The input matrix C is changed to Q' * C if desired. * * See "Computing Small Singular Values of Bidiagonal Matrices With * Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan, * LAPACK Working Note #3, for a detailed description of the algorithm. * * Arguments * ========= * * UPLO (input) CHARACTER*1 * On entry, UPLO specifies whether the input bidiagonal matrix * is upper or lower bidiagonal, and wether it is square are * not. * UPLO = 'U' or 'u' B is upper bidiagonal. * UPLO = 'L' or 'l' B is lower bidiagonal. * * SQRE (input) INTEGER * = 0: then the input matrix is N-by-N. * = 1: then the input matrix is N-by-(N+1) if UPLU = 'U' and * (N+1)-by-N if UPLU = 'L'. * * The bidiagonal matrix has * N = NL + NR + 1 rows and * M = N + SQRE >= N columns. * * N (input) INTEGER * On entry, N specifies the number of rows and columns * in the matrix. N must be at least 0. * * NCVT (input) INTEGER * On entry, NCVT specifies the number of columns of * the matrix VT. NCVT must be at least 0. * * NRU (input) INTEGER * On entry, NRU specifies the number of rows of * the matrix U. NRU must be at least 0. * * NCC (input) INTEGER * On entry, NCC specifies the number of columns of * the matrix C. NCC must be at least 0. * * D (input/output) DOUBLE PRECISION array, dimension (N) * On entry, D contains the diagonal entries of the * bidiagonal matrix whose SVD is desired. On normal exit, * D contains the singular values in ascending order. * * E (input/output) DOUBLE PRECISION array. * dimension is (N-1) if SQRE = 0 and N if SQRE = 1. * On entry, the entries of E contain the offdiagonal entries * of the bidiagonal matrix whose SVD is desired. On normal * exit, E will contain 0. If the algorithm does not converge, * D and E will contain the diagonal and superdiagonal entries * of a bidiagonal matrix orthogonally equivalent to the one * given as input. * * VT (input/output) DOUBLE PRECISION array, dimension (LDVT, NCVT) * On entry, contains a matrix which on exit has been * premultiplied by P', dimension N-by-NCVT if SQRE = 0 * and (N+1)-by-NCVT if SQRE = 1 (not referenced if NCVT=0). * * LDVT (input) INTEGER * On entry, LDVT specifies the leading dimension of VT as * declared in the calling (sub) program. LDVT must be at * least 1. If NCVT is nonzero LDVT must also be at least N. * * U (input/output) DOUBLE PRECISION array, dimension (LDU, N) * On entry, contains a matrix which on exit has been * postmultiplied by Q, dimension NRU-by-N if SQRE = 0 * and NRU-by-(N+1) if SQRE = 1 (not referenced if NRU=0). * * LDU (input) INTEGER * On entry, LDU specifies the leading dimension of U as * declared in the calling (sub) program. LDU must be at * least max( 1, NRU ) . * * C (input/output) DOUBLE PRECISION array, dimension (LDC, NCC) * On entry, contains an N-by-NCC matrix which on exit * has been premultiplied by Q' dimension N-by-NCC if SQRE = 0 * and (N+1)-by-NCC if SQRE = 1 (not referenced if NCC=0). * * LDC (input) INTEGER * On entry, LDC specifies the leading dimension of C as * declared in the calling (sub) program. LDC must be at * least 1. If NCC is nonzero, LDC must also be at least N. * * WORK (workspace) DOUBLE PRECISION array, dimension (4*N) * Workspace. Only referenced if one of NCVT, NRU, or NCC is * nonzero, and if N is at least 2. * * INFO (output) INTEGER * On exit, a value of 0 indicates a successful exit. * If INFO < 0, argument number -INFO is illegal. * If INFO > 0, the algorithm did not converge, and INFO * specifies how many superdiagonals did not converge. * * Further Details * =============== * * Based on contributions by * Ming Gu and Huan Ren, Computer Science Division, University of * California at Berkeley, USA * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) * .. * .. Local Scalars .. LOGICAL ROTATE INTEGER I, ISUB, IUPLO, J, NP1, SQRE1 DOUBLE PRECISION CS, R, SMIN, SN * .. * .. External Subroutines .. EXTERNAL DBDSQR, DLARTG, DLASR, DSWAP, XERBLA * .. * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IUPLO = 0 IF( LSAME( UPLO, 'U' ) ) \$ IUPLO = 1 IF( LSAME( UPLO, 'L' ) ) \$ IUPLO = 2 IF( IUPLO.EQ.0 ) THEN INFO = -1 ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN INFO = -2 ELSE IF( N.LT.0 ) THEN INFO = -3 ELSE IF( NCVT.LT.0 ) THEN INFO = -4 ELSE IF( NRU.LT.0 ) THEN INFO = -5 ELSE IF( NCC.LT.0 ) THEN INFO = -6 ELSE IF( ( NCVT.EQ.0 .AND. LDVT.LT.1 ) .OR. \$ ( NCVT.GT.0 .AND. LDVT.LT.MAX( 1, N ) ) ) THEN INFO = -10 ELSE IF( LDU.LT.MAX( 1, NRU ) ) THEN INFO = -12 ELSE IF( ( NCC.EQ.0 .AND. LDC.LT.1 ) .OR. \$ ( NCC.GT.0 .AND. LDC.LT.MAX( 1, N ) ) ) THEN INFO = -14 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DLASDQ', -INFO ) RETURN END IF IF( N.EQ.0 ) \$ RETURN * * ROTATE is true if any singular vectors desired, false otherwise * ROTATE = ( NCVT.GT.0 ) .OR. ( NRU.GT.0 ) .OR. ( NCC.GT.0 ) NP1 = N + 1 SQRE1 = SQRE * * If matrix non-square upper bidiagonal, rotate to be lower * bidiagonal. The rotations are on the right. * IF( ( IUPLO.EQ.1 ) .AND. ( SQRE1.EQ.1 ) ) 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 ) IF( ROTATE ) THEN WORK( I ) = CS WORK( N+I ) = SN END IF 10 CONTINUE CALL DLARTG( D( N ), E( N ), CS, SN, R ) D( N ) = R E( N ) = ZERO IF( ROTATE ) THEN WORK( N ) = CS WORK( N+N ) = SN END IF IUPLO = 2 SQRE1 = 0 * * Update singular vectors if desired. * IF( NCVT.GT.0 ) \$ CALL DLASR( 'L', 'V', 'F', NP1, NCVT, WORK( 1 ), \$ WORK( NP1 ), VT, LDVT ) END IF * * If matrix lower bidiagonal, rotate to be upper bidiagonal * by applying Givens rotations on the left. * IF( IUPLO.EQ.2 ) THEN DO 20 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 ) IF( ROTATE ) THEN WORK( I ) = CS WORK( N+I ) = SN END IF 20 CONTINUE * * If matrix (N+1)-by-N lower bidiagonal, one additional * rotation is needed. * IF( SQRE1.EQ.1 ) THEN CALL DLARTG( D( N ), E( N ), CS, SN, R ) D( N ) = R IF( ROTATE ) THEN WORK( N ) = CS WORK( N+N ) = SN END IF END IF * * Update singular vectors if desired. * IF( NRU.GT.0 ) THEN IF( SQRE1.EQ.0 ) THEN CALL DLASR( 'R', 'V', 'F', NRU, N, WORK( 1 ), \$ WORK( NP1 ), U, LDU ) ELSE CALL DLASR( 'R', 'V', 'F', NRU, NP1, WORK( 1 ), \$ WORK( NP1 ), U, LDU ) END IF END IF IF( NCC.GT.0 ) THEN IF( SQRE1.EQ.0 ) THEN CALL DLASR( 'L', 'V', 'F', N, NCC, WORK( 1 ), \$ WORK( NP1 ), C, LDC ) ELSE CALL DLASR( 'L', 'V', 'F', NP1, NCC, WORK( 1 ), \$ WORK( NP1 ), C, LDC ) END IF END IF END IF * * Call DBDSQR to compute the SVD of the reduced real * N-by-N upper bidiagonal matrix. * CALL DBDSQR( 'U', N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, \$ LDC, WORK, INFO ) * * Sort the singular values into ascending order (insertion sort on * singular values, but only one transposition per singular vector) * DO 40 I = 1, N * * Scan for smallest D(I). * ISUB = I SMIN = D( I ) DO 30 J = I + 1, N IF( D( J ).LT.SMIN ) THEN ISUB = J SMIN = D( J ) END IF 30 CONTINUE IF( ISUB.NE.I ) THEN * * Swap singular values and vectors. * D( ISUB ) = D( I ) D( I ) = SMIN IF( NCVT.GT.0 ) \$ CALL DSWAP( NCVT, VT( ISUB, 1 ), LDVT, VT( I, 1 ), LDVT ) IF( NRU.GT.0 ) \$ CALL DSWAP( NRU, U( 1, ISUB ), 1, U( 1, I ), 1 ) IF( NCC.GT.0 ) \$ CALL DSWAP( NCC, C( ISUB, 1 ), LDC, C( I, 1 ), LDC ) END IF 40 CONTINUE * RETURN * * End of DLASDQ * END