![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DGGSVD( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, 2: $ LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, 3: $ IWORK, INFO ) 4: * 5: * -- LAPACK driver routine (version 3.2) -- 6: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 7: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 8: * November 2006 9: * 10: * .. Scalar Arguments .. 11: CHARACTER JOBQ, JOBU, JOBV 12: INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P 13: * .. 14: * .. Array Arguments .. 15: INTEGER IWORK( * ) 16: DOUBLE PRECISION A( LDA, * ), ALPHA( * ), B( LDB, * ), 17: $ BETA( * ), Q( LDQ, * ), U( LDU, * ), 18: $ V( LDV, * ), WORK( * ) 19: * .. 20: * 21: * Purpose 22: * ======= 23: * 24: * DGGSVD computes the generalized singular value decomposition (GSVD) 25: * of an M-by-N real matrix A and P-by-N real matrix B: 26: * 27: * U'*A*Q = D1*( 0 R ), V'*B*Q = D2*( 0 R ) 28: * 29: * where U, V and Q are orthogonal matrices, and Z' is the transpose 30: * of Z. Let K+L = the effective numerical rank of the matrix (A',B')', 31: * then R is a K+L-by-K+L nonsingular upper triangular matrix, D1 and 32: * D2 are M-by-(K+L) and P-by-(K+L) "diagonal" matrices and of the 33: * following structures, respectively: 34: * 35: * If M-K-L >= 0, 36: * 37: * K L 38: * D1 = K ( I 0 ) 39: * L ( 0 C ) 40: * M-K-L ( 0 0 ) 41: * 42: * K L 43: * D2 = L ( 0 S ) 44: * P-L ( 0 0 ) 45: * 46: * N-K-L K L 47: * ( 0 R ) = K ( 0 R11 R12 ) 48: * L ( 0 0 R22 ) 49: * 50: * where 51: * 52: * C = diag( ALPHA(K+1), ... , ALPHA(K+L) ), 53: * S = diag( BETA(K+1), ... , BETA(K+L) ), 54: * C**2 + S**2 = I. 55: * 56: * R is stored in A(1:K+L,N-K-L+1:N) on exit. 57: * 58: * If M-K-L < 0, 59: * 60: * K M-K K+L-M 61: * D1 = K ( I 0 0 ) 62: * M-K ( 0 C 0 ) 63: * 64: * K M-K K+L-M 65: * D2 = M-K ( 0 S 0 ) 66: * K+L-M ( 0 0 I ) 67: * P-L ( 0 0 0 ) 68: * 69: * N-K-L K M-K K+L-M 70: * ( 0 R ) = K ( 0 R11 R12 R13 ) 71: * M-K ( 0 0 R22 R23 ) 72: * K+L-M ( 0 0 0 R33 ) 73: * 74: * where 75: * 76: * C = diag( ALPHA(K+1), ... , ALPHA(M) ), 77: * S = diag( BETA(K+1), ... , BETA(M) ), 78: * C**2 + S**2 = I. 79: * 80: * (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), and R33 is stored 81: * ( 0 R22 R23 ) 82: * in B(M-K+1:L,N+M-K-L+1:N) on exit. 83: * 84: * The routine computes C, S, R, and optionally the orthogonal 85: * transformation matrices U, V and Q. 86: * 87: * In particular, if B is an N-by-N nonsingular matrix, then the GSVD of 88: * A and B implicitly gives the SVD of A*inv(B): 89: * A*inv(B) = U*(D1*inv(D2))*V'. 90: * If ( A',B')' has orthonormal columns, then the GSVD of A and B is 91: * also equal to the CS decomposition of A and B. Furthermore, the GSVD 92: * can be used to derive the solution of the eigenvalue problem: 93: * A'*A x = lambda* B'*B x. 94: * In some literature, the GSVD of A and B is presented in the form 95: * U'*A*X = ( 0 D1 ), V'*B*X = ( 0 D2 ) 96: * where U and V are orthogonal and X is nonsingular, D1 and D2 are 97: * ``diagonal''. The former GSVD form can be converted to the latter 98: * form by taking the nonsingular matrix X as 99: * 100: * X = Q*( I 0 ) 101: * ( 0 inv(R) ). 102: * 103: * Arguments 104: * ========= 105: * 106: * JOBU (input) CHARACTER*1 107: * = 'U': Orthogonal matrix U is computed; 108: * = 'N': U is not computed. 109: * 110: * JOBV (input) CHARACTER*1 111: * = 'V': Orthogonal matrix V is computed; 112: * = 'N': V is not computed. 113: * 114: * JOBQ (input) CHARACTER*1 115: * = 'Q': Orthogonal matrix Q is computed; 116: * = 'N': Q is not computed. 117: * 118: * M (input) INTEGER 119: * The number of rows of the matrix A. M >= 0. 120: * 121: * N (input) INTEGER 122: * The number of columns of the matrices A and B. N >= 0. 123: * 124: * P (input) INTEGER 125: * The number of rows of the matrix B. P >= 0. 126: * 127: * K (output) INTEGER 128: * L (output) INTEGER 129: * On exit, K and L specify the dimension of the subblocks 130: * described in the Purpose section. 131: * K + L = effective numerical rank of (A',B')'. 132: * 133: * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 134: * On entry, the M-by-N matrix A. 135: * On exit, A contains the triangular matrix R, or part of R. 136: * See Purpose for details. 137: * 138: * LDA (input) INTEGER 139: * The leading dimension of the array A. LDA >= max(1,M). 140: * 141: * B (input/output) DOUBLE PRECISION array, dimension (LDB,N) 142: * On entry, the P-by-N matrix B. 143: * On exit, B contains the triangular matrix R if M-K-L < 0. 144: * See Purpose for details. 145: * 146: * LDB (input) INTEGER 147: * The leading dimension of the array B. LDB >= max(1,P). 148: * 149: * ALPHA (output) DOUBLE PRECISION array, dimension (N) 150: * BETA (output) DOUBLE PRECISION array, dimension (N) 151: * On exit, ALPHA and BETA contain the generalized singular 152: * value pairs of A and B; 153: * ALPHA(1:K) = 1, 154: * BETA(1:K) = 0, 155: * and if M-K-L >= 0, 156: * ALPHA(K+1:K+L) = C, 157: * BETA(K+1:K+L) = S, 158: * or if M-K-L < 0, 159: * ALPHA(K+1:M)=C, ALPHA(M+1:K+L)=0 160: * BETA(K+1:M) =S, BETA(M+1:K+L) =1 161: * and 162: * ALPHA(K+L+1:N) = 0 163: * BETA(K+L+1:N) = 0 164: * 165: * U (output) DOUBLE PRECISION array, dimension (LDU,M) 166: * If JOBU = 'U', U contains the M-by-M orthogonal matrix U. 167: * If JOBU = 'N', U is not referenced. 168: * 169: * LDU (input) INTEGER 170: * The leading dimension of the array U. LDU >= max(1,M) if 171: * JOBU = 'U'; LDU >= 1 otherwise. 172: * 173: * V (output) DOUBLE PRECISION array, dimension (LDV,P) 174: * If JOBV = 'V', V contains the P-by-P orthogonal matrix V. 175: * If JOBV = 'N', V is not referenced. 176: * 177: * LDV (input) INTEGER 178: * The leading dimension of the array V. LDV >= max(1,P) if 179: * JOBV = 'V'; LDV >= 1 otherwise. 180: * 181: * Q (output) DOUBLE PRECISION array, dimension (LDQ,N) 182: * If JOBQ = 'Q', Q contains the N-by-N orthogonal matrix Q. 183: * If JOBQ = 'N', Q is not referenced. 184: * 185: * LDQ (input) INTEGER 186: * The leading dimension of the array Q. LDQ >= max(1,N) if 187: * JOBQ = 'Q'; LDQ >= 1 otherwise. 188: * 189: * WORK (workspace) DOUBLE PRECISION array, 190: * dimension (max(3*N,M,P)+N) 191: * 192: * IWORK (workspace/output) INTEGER array, dimension (N) 193: * On exit, IWORK stores the sorting information. More 194: * precisely, the following loop will sort ALPHA 195: * for I = K+1, min(M,K+L) 196: * swap ALPHA(I) and ALPHA(IWORK(I)) 197: * endfor 198: * such that ALPHA(1) >= ALPHA(2) >= ... >= ALPHA(N). 199: * 200: * INFO (output) INTEGER 201: * = 0: successful exit 202: * < 0: if INFO = -i, the i-th argument had an illegal value. 203: * > 0: if INFO = 1, the Jacobi-type procedure failed to 204: * converge. For further details, see subroutine DTGSJA. 205: * 206: * Internal Parameters 207: * =================== 208: * 209: * TOLA DOUBLE PRECISION 210: * TOLB DOUBLE PRECISION 211: * TOLA and TOLB are the thresholds to determine the effective 212: * rank of (A',B')'. Generally, they are set to 213: * TOLA = MAX(M,N)*norm(A)*MAZHEPS, 214: * TOLB = MAX(P,N)*norm(B)*MAZHEPS. 215: * The size of TOLA and TOLB may affect the size of backward 216: * errors of the decomposition. 217: * 218: * Further Details 219: * =============== 220: * 221: * 2-96 Based on modifications by 222: * Ming Gu and Huan Ren, Computer Science Division, University of 223: * California at Berkeley, USA 224: * 225: * ===================================================================== 226: * 227: * .. Local Scalars .. 228: LOGICAL WANTQ, WANTU, WANTV 229: INTEGER I, IBND, ISUB, J, NCYCLE 230: DOUBLE PRECISION ANORM, BNORM, SMAX, TEMP, TOLA, TOLB, ULP, UNFL 231: * .. 232: * .. External Functions .. 233: LOGICAL LSAME 234: DOUBLE PRECISION DLAMCH, DLANGE 235: EXTERNAL LSAME, DLAMCH, DLANGE 236: * .. 237: * .. External Subroutines .. 238: EXTERNAL DCOPY, DGGSVP, DTGSJA, XERBLA 239: * .. 240: * .. Intrinsic Functions .. 241: INTRINSIC MAX, MIN 242: * .. 243: * .. Executable Statements .. 244: * 245: * Test the input parameters 246: * 247: WANTU = LSAME( JOBU, 'U' ) 248: WANTV = LSAME( JOBV, 'V' ) 249: WANTQ = LSAME( JOBQ, 'Q' ) 250: * 251: INFO = 0 252: IF( .NOT.( WANTU .OR. LSAME( JOBU, 'N' ) ) ) THEN 253: INFO = -1 254: ELSE IF( .NOT.( WANTV .OR. LSAME( JOBV, 'N' ) ) ) THEN 255: INFO = -2 256: ELSE IF( .NOT.( WANTQ .OR. LSAME( JOBQ, 'N' ) ) ) THEN 257: INFO = -3 258: ELSE IF( M.LT.0 ) THEN 259: INFO = -4 260: ELSE IF( N.LT.0 ) THEN 261: INFO = -5 262: ELSE IF( P.LT.0 ) THEN 263: INFO = -6 264: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 265: INFO = -10 266: ELSE IF( LDB.LT.MAX( 1, P ) ) THEN 267: INFO = -12 268: ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN 269: INFO = -16 270: ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN 271: INFO = -18 272: ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN 273: INFO = -20 274: END IF 275: IF( INFO.NE.0 ) THEN 276: CALL XERBLA( 'DGGSVD', -INFO ) 277: RETURN 278: END IF 279: * 280: * Compute the Frobenius norm of matrices A and B 281: * 282: ANORM = DLANGE( '1', M, N, A, LDA, WORK ) 283: BNORM = DLANGE( '1', P, N, B, LDB, WORK ) 284: * 285: * Get machine precision and set up threshold for determining 286: * the effective numerical rank of the matrices A and B. 287: * 288: ULP = DLAMCH( 'Precision' ) 289: UNFL = DLAMCH( 'Safe Minimum' ) 290: TOLA = MAX( M, N )*MAX( ANORM, UNFL )*ULP 291: TOLB = MAX( P, N )*MAX( BNORM, UNFL )*ULP 292: * 293: * Preprocessing 294: * 295: CALL DGGSVP( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, TOLA, 296: $ TOLB, K, L, U, LDU, V, LDV, Q, LDQ, IWORK, WORK, 297: $ WORK( N+1 ), INFO ) 298: * 299: * Compute the GSVD of two upper "triangular" matrices 300: * 301: CALL DTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, LDB, 302: $ TOLA, TOLB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, 303: $ WORK, NCYCLE, INFO ) 304: * 305: * Sort the singular values and store the pivot indices in IWORK 306: * Copy ALPHA to WORK, then sort ALPHA in WORK 307: * 308: CALL DCOPY( N, ALPHA, 1, WORK, 1 ) 309: IBND = MIN( L, M-K ) 310: DO 20 I = 1, IBND 311: * 312: * Scan for largest ALPHA(K+I) 313: * 314: ISUB = I 315: SMAX = WORK( K+I ) 316: DO 10 J = I + 1, IBND 317: TEMP = WORK( K+J ) 318: IF( TEMP.GT.SMAX ) THEN 319: ISUB = J 320: SMAX = TEMP 321: END IF 322: 10 CONTINUE 323: IF( ISUB.NE.I ) THEN 324: WORK( K+ISUB ) = WORK( K+I ) 325: WORK( K+I ) = SMAX 326: IWORK( K+I ) = K + ISUB 327: ELSE 328: IWORK( K+I ) = K + I 329: END IF 330: 20 CONTINUE 331: * 332: RETURN 333: * 334: * End of DGGSVD 335: * 336: END