Annotation of rpl/lapack/lapack/dlasdq.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DLASDQ( UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT,
! 2: $ U, LDU, C, LDC, WORK, INFO )
! 3: *
! 4: * -- LAPACK auxiliary routine (version 3.2) --
! 5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 7: * November 2006
! 8: *
! 9: * .. Scalar Arguments ..
! 10: CHARACTER UPLO
! 11: INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU, SQRE
! 12: * ..
! 13: * .. Array Arguments ..
! 14: DOUBLE PRECISION C( LDC, * ), D( * ), E( * ), U( LDU, * ),
! 15: $ VT( LDVT, * ), WORK( * )
! 16: * ..
! 17: *
! 18: * Purpose
! 19: * =======
! 20: *
! 21: * DLASDQ computes the singular value decomposition (SVD) of a real
! 22: * (upper or lower) bidiagonal matrix with diagonal D and offdiagonal
! 23: * E, accumulating the transformations if desired. Letting B denote
! 24: * the input bidiagonal matrix, the algorithm computes orthogonal
! 25: * matrices Q and P such that B = Q * S * P' (P' denotes the transpose
! 26: * of P). The singular values S are overwritten on D.
! 27: *
! 28: * The input matrix U is changed to U * Q if desired.
! 29: * The input matrix VT is changed to P' * VT if desired.
! 30: * The input matrix C is changed to Q' * C if desired.
! 31: *
! 32: * See "Computing Small Singular Values of Bidiagonal Matrices With
! 33: * Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
! 34: * LAPACK Working Note #3, for a detailed description of the algorithm.
! 35: *
! 36: * Arguments
! 37: * =========
! 38: *
! 39: * UPLO (input) CHARACTER*1
! 40: * On entry, UPLO specifies whether the input bidiagonal matrix
! 41: * is upper or lower bidiagonal, and wether it is square are
! 42: * not.
! 43: * UPLO = 'U' or 'u' B is upper bidiagonal.
! 44: * UPLO = 'L' or 'l' B is lower bidiagonal.
! 45: *
! 46: * SQRE (input) INTEGER
! 47: * = 0: then the input matrix is N-by-N.
! 48: * = 1: then the input matrix is N-by-(N+1) if UPLU = 'U' and
! 49: * (N+1)-by-N if UPLU = 'L'.
! 50: *
! 51: * The bidiagonal matrix has
! 52: * N = NL + NR + 1 rows and
! 53: * M = N + SQRE >= N columns.
! 54: *
! 55: * N (input) INTEGER
! 56: * On entry, N specifies the number of rows and columns
! 57: * in the matrix. N must be at least 0.
! 58: *
! 59: * NCVT (input) INTEGER
! 60: * On entry, NCVT specifies the number of columns of
! 61: * the matrix VT. NCVT must be at least 0.
! 62: *
! 63: * NRU (input) INTEGER
! 64: * On entry, NRU specifies the number of rows of
! 65: * the matrix U. NRU must be at least 0.
! 66: *
! 67: * NCC (input) INTEGER
! 68: * On entry, NCC specifies the number of columns of
! 69: * the matrix C. NCC must be at least 0.
! 70: *
! 71: * D (input/output) DOUBLE PRECISION array, dimension (N)
! 72: * On entry, D contains the diagonal entries of the
! 73: * bidiagonal matrix whose SVD is desired. On normal exit,
! 74: * D contains the singular values in ascending order.
! 75: *
! 76: * E (input/output) DOUBLE PRECISION array.
! 77: * dimension is (N-1) if SQRE = 0 and N if SQRE = 1.
! 78: * On entry, the entries of E contain the offdiagonal entries
! 79: * of the bidiagonal matrix whose SVD is desired. On normal
! 80: * exit, E will contain 0. If the algorithm does not converge,
! 81: * D and E will contain the diagonal and superdiagonal entries
! 82: * of a bidiagonal matrix orthogonally equivalent to the one
! 83: * given as input.
! 84: *
! 85: * VT (input/output) DOUBLE PRECISION array, dimension (LDVT, NCVT)
! 86: * On entry, contains a matrix which on exit has been
! 87: * premultiplied by P', dimension N-by-NCVT if SQRE = 0
! 88: * and (N+1)-by-NCVT if SQRE = 1 (not referenced if NCVT=0).
! 89: *
! 90: * LDVT (input) INTEGER
! 91: * On entry, LDVT specifies the leading dimension of VT as
! 92: * declared in the calling (sub) program. LDVT must be at
! 93: * least 1. If NCVT is nonzero LDVT must also be at least N.
! 94: *
! 95: * U (input/output) DOUBLE PRECISION array, dimension (LDU, N)
! 96: * On entry, contains a matrix which on exit has been
! 97: * postmultiplied by Q, dimension NRU-by-N if SQRE = 0
! 98: * and NRU-by-(N+1) if SQRE = 1 (not referenced if NRU=0).
! 99: *
! 100: * LDU (input) INTEGER
! 101: * On entry, LDU specifies the leading dimension of U as
! 102: * declared in the calling (sub) program. LDU must be at
! 103: * least max( 1, NRU ) .
! 104: *
! 105: * C (input/output) DOUBLE PRECISION array, dimension (LDC, NCC)
! 106: * On entry, contains an N-by-NCC matrix which on exit
! 107: * has been premultiplied by Q' dimension N-by-NCC if SQRE = 0
! 108: * and (N+1)-by-NCC if SQRE = 1 (not referenced if NCC=0).
! 109: *
! 110: * LDC (input) INTEGER
! 111: * On entry, LDC specifies the leading dimension of C as
! 112: * declared in the calling (sub) program. LDC must be at
! 113: * least 1. If NCC is nonzero, LDC must also be at least N.
! 114: *
! 115: * WORK (workspace) DOUBLE PRECISION array, dimension (4*N)
! 116: * Workspace. Only referenced if one of NCVT, NRU, or NCC is
! 117: * nonzero, and if N is at least 2.
! 118: *
! 119: * INFO (output) INTEGER
! 120: * On exit, a value of 0 indicates a successful exit.
! 121: * If INFO < 0, argument number -INFO is illegal.
! 122: * If INFO > 0, the algorithm did not converge, and INFO
! 123: * specifies how many superdiagonals did not converge.
! 124: *
! 125: * Further Details
! 126: * ===============
! 127: *
! 128: * Based on contributions by
! 129: * Ming Gu and Huan Ren, Computer Science Division, University of
! 130: * California at Berkeley, USA
! 131: *
! 132: * =====================================================================
! 133: *
! 134: * .. Parameters ..
! 135: DOUBLE PRECISION ZERO
! 136: PARAMETER ( ZERO = 0.0D+0 )
! 137: * ..
! 138: * .. Local Scalars ..
! 139: LOGICAL ROTATE
! 140: INTEGER I, ISUB, IUPLO, J, NP1, SQRE1
! 141: DOUBLE PRECISION CS, R, SMIN, SN
! 142: * ..
! 143: * .. External Subroutines ..
! 144: EXTERNAL DBDSQR, DLARTG, DLASR, DSWAP, XERBLA
! 145: * ..
! 146: * .. External Functions ..
! 147: LOGICAL LSAME
! 148: EXTERNAL LSAME
! 149: * ..
! 150: * .. Intrinsic Functions ..
! 151: INTRINSIC MAX
! 152: * ..
! 153: * .. Executable Statements ..
! 154: *
! 155: * Test the input parameters.
! 156: *
! 157: INFO = 0
! 158: IUPLO = 0
! 159: IF( LSAME( UPLO, 'U' ) )
! 160: $ IUPLO = 1
! 161: IF( LSAME( UPLO, 'L' ) )
! 162: $ IUPLO = 2
! 163: IF( IUPLO.EQ.0 ) THEN
! 164: INFO = -1
! 165: ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
! 166: INFO = -2
! 167: ELSE IF( N.LT.0 ) THEN
! 168: INFO = -3
! 169: ELSE IF( NCVT.LT.0 ) THEN
! 170: INFO = -4
! 171: ELSE IF( NRU.LT.0 ) THEN
! 172: INFO = -5
! 173: ELSE IF( NCC.LT.0 ) THEN
! 174: INFO = -6
! 175: ELSE IF( ( NCVT.EQ.0 .AND. LDVT.LT.1 ) .OR.
! 176: $ ( NCVT.GT.0 .AND. LDVT.LT.MAX( 1, N ) ) ) THEN
! 177: INFO = -10
! 178: ELSE IF( LDU.LT.MAX( 1, NRU ) ) THEN
! 179: INFO = -12
! 180: ELSE IF( ( NCC.EQ.0 .AND. LDC.LT.1 ) .OR.
! 181: $ ( NCC.GT.0 .AND. LDC.LT.MAX( 1, N ) ) ) THEN
! 182: INFO = -14
! 183: END IF
! 184: IF( INFO.NE.0 ) THEN
! 185: CALL XERBLA( 'DLASDQ', -INFO )
! 186: RETURN
! 187: END IF
! 188: IF( N.EQ.0 )
! 189: $ RETURN
! 190: *
! 191: * ROTATE is true if any singular vectors desired, false otherwise
! 192: *
! 193: ROTATE = ( NCVT.GT.0 ) .OR. ( NRU.GT.0 ) .OR. ( NCC.GT.0 )
! 194: NP1 = N + 1
! 195: SQRE1 = SQRE
! 196: *
! 197: * If matrix non-square upper bidiagonal, rotate to be lower
! 198: * bidiagonal. The rotations are on the right.
! 199: *
! 200: IF( ( IUPLO.EQ.1 ) .AND. ( SQRE1.EQ.1 ) ) THEN
! 201: DO 10 I = 1, N - 1
! 202: CALL DLARTG( D( I ), E( I ), CS, SN, R )
! 203: D( I ) = R
! 204: E( I ) = SN*D( I+1 )
! 205: D( I+1 ) = CS*D( I+1 )
! 206: IF( ROTATE ) THEN
! 207: WORK( I ) = CS
! 208: WORK( N+I ) = SN
! 209: END IF
! 210: 10 CONTINUE
! 211: CALL DLARTG( D( N ), E( N ), CS, SN, R )
! 212: D( N ) = R
! 213: E( N ) = ZERO
! 214: IF( ROTATE ) THEN
! 215: WORK( N ) = CS
! 216: WORK( N+N ) = SN
! 217: END IF
! 218: IUPLO = 2
! 219: SQRE1 = 0
! 220: *
! 221: * Update singular vectors if desired.
! 222: *
! 223: IF( NCVT.GT.0 )
! 224: $ CALL DLASR( 'L', 'V', 'F', NP1, NCVT, WORK( 1 ),
! 225: $ WORK( NP1 ), VT, LDVT )
! 226: END IF
! 227: *
! 228: * If matrix lower bidiagonal, rotate to be upper bidiagonal
! 229: * by applying Givens rotations on the left.
! 230: *
! 231: IF( IUPLO.EQ.2 ) THEN
! 232: DO 20 I = 1, N - 1
! 233: CALL DLARTG( D( I ), E( I ), CS, SN, R )
! 234: D( I ) = R
! 235: E( I ) = SN*D( I+1 )
! 236: D( I+1 ) = CS*D( I+1 )
! 237: IF( ROTATE ) THEN
! 238: WORK( I ) = CS
! 239: WORK( N+I ) = SN
! 240: END IF
! 241: 20 CONTINUE
! 242: *
! 243: * If matrix (N+1)-by-N lower bidiagonal, one additional
! 244: * rotation is needed.
! 245: *
! 246: IF( SQRE1.EQ.1 ) THEN
! 247: CALL DLARTG( D( N ), E( N ), CS, SN, R )
! 248: D( N ) = R
! 249: IF( ROTATE ) THEN
! 250: WORK( N ) = CS
! 251: WORK( N+N ) = SN
! 252: END IF
! 253: END IF
! 254: *
! 255: * Update singular vectors if desired.
! 256: *
! 257: IF( NRU.GT.0 ) THEN
! 258: IF( SQRE1.EQ.0 ) THEN
! 259: CALL DLASR( 'R', 'V', 'F', NRU, N, WORK( 1 ),
! 260: $ WORK( NP1 ), U, LDU )
! 261: ELSE
! 262: CALL DLASR( 'R', 'V', 'F', NRU, NP1, WORK( 1 ),
! 263: $ WORK( NP1 ), U, LDU )
! 264: END IF
! 265: END IF
! 266: IF( NCC.GT.0 ) THEN
! 267: IF( SQRE1.EQ.0 ) THEN
! 268: CALL DLASR( 'L', 'V', 'F', N, NCC, WORK( 1 ),
! 269: $ WORK( NP1 ), C, LDC )
! 270: ELSE
! 271: CALL DLASR( 'L', 'V', 'F', NP1, NCC, WORK( 1 ),
! 272: $ WORK( NP1 ), C, LDC )
! 273: END IF
! 274: END IF
! 275: END IF
! 276: *
! 277: * Call DBDSQR to compute the SVD of the reduced real
! 278: * N-by-N upper bidiagonal matrix.
! 279: *
! 280: CALL DBDSQR( 'U', N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C,
! 281: $ LDC, WORK, INFO )
! 282: *
! 283: * Sort the singular values into ascending order (insertion sort on
! 284: * singular values, but only one transposition per singular vector)
! 285: *
! 286: DO 40 I = 1, N
! 287: *
! 288: * Scan for smallest D(I).
! 289: *
! 290: ISUB = I
! 291: SMIN = D( I )
! 292: DO 30 J = I + 1, N
! 293: IF( D( J ).LT.SMIN ) THEN
! 294: ISUB = J
! 295: SMIN = D( J )
! 296: END IF
! 297: 30 CONTINUE
! 298: IF( ISUB.NE.I ) THEN
! 299: *
! 300: * Swap singular values and vectors.
! 301: *
! 302: D( ISUB ) = D( I )
! 303: D( I ) = SMIN
! 304: IF( NCVT.GT.0 )
! 305: $ CALL DSWAP( NCVT, VT( ISUB, 1 ), LDVT, VT( I, 1 ), LDVT )
! 306: IF( NRU.GT.0 )
! 307: $ CALL DSWAP( NRU, U( 1, ISUB ), 1, U( 1, I ), 1 )
! 308: IF( NCC.GT.0 )
! 309: $ CALL DSWAP( NCC, C( ISUB, 1 ), LDC, C( I, 1 ), LDC )
! 310: END IF
! 311: 40 CONTINUE
! 312: *
! 313: RETURN
! 314: *
! 315: * End of DLASDQ
! 316: *
! 317: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasdq.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack