![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
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