![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZGGSVD( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, 2: $ LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, 3: $ RWORK, 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 ALPHA( * ), BETA( * ), RWORK( * ) 17: COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ), 18: $ U( LDU, * ), V( LDV, * ), WORK( * ) 19: * .. 20: * 21: * Purpose 22: * ======= 23: * 24: * ZGGSVD computes the generalized singular value decomposition (GSVD) 25: * of an M-by-N complex matrix A and P-by-N complex matrix B: 26: * 27: * U'*A*Q = D1*( 0 R ), V'*B*Q = D2*( 0 R ) 28: * 29: * where U, V and Q are unitary matrices, and Z' means the conjugate 30: * transpose of Z. Let K+L = the effective numerical rank of the 31: * matrix (A',B')', then R is a (K+L)-by-(K+L) nonsingular upper 32: * triangular matrix, D1 and D2 are M-by-(K+L) and P-by-(K+L) "diagonal" 33: * matrices and of the 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: * where 50: * 51: * C = diag( ALPHA(K+1), ... , ALPHA(K+L) ), 52: * S = diag( BETA(K+1), ... , BETA(K+L) ), 53: * C**2 + S**2 = I. 54: * 55: * R is stored in A(1:K+L,N-K-L+1:N) on exit. 56: * 57: * If M-K-L < 0, 58: * 59: * K M-K K+L-M 60: * D1 = K ( I 0 0 ) 61: * M-K ( 0 C 0 ) 62: * 63: * K M-K K+L-M 64: * D2 = M-K ( 0 S 0 ) 65: * K+L-M ( 0 0 I ) 66: * P-L ( 0 0 0 ) 67: * 68: * N-K-L K M-K K+L-M 69: * ( 0 R ) = K ( 0 R11 R12 R13 ) 70: * M-K ( 0 0 R22 R23 ) 71: * K+L-M ( 0 0 0 R33 ) 72: * 73: * where 74: * 75: * C = diag( ALPHA(K+1), ... , ALPHA(M) ), 76: * S = diag( BETA(K+1), ... , BETA(M) ), 77: * C**2 + S**2 = I. 78: * 79: * (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), and R33 is stored 80: * ( 0 R22 R23 ) 81: * in B(M-K+1:L,N+M-K-L+1:N) on exit. 82: * 83: * The routine computes C, S, R, and optionally the unitary 84: * transformation matrices U, V and Q. 85: * 86: * In particular, if B is an N-by-N nonsingular matrix, then the GSVD of 87: * A and B implicitly gives the SVD of A*inv(B): 88: * A*inv(B) = U*(D1*inv(D2))*V'. 89: * If ( A',B')' has orthnormal columns, then the GSVD of A and B is also 90: * equal to the CS decomposition of A and B. Furthermore, the GSVD can 91: * be used to derive the solution of the eigenvalue problem: 92: * A'*A x = lambda* B'*B x. 93: * In some literature, the GSVD of A and B is presented in the form 94: * U'*A*X = ( 0 D1 ), V'*B*X = ( 0 D2 ) 95: * where U and V are orthogonal and X is nonsingular, and D1 and D2 are 96: * ``diagonal''. The former GSVD form can be converted to the latter 97: * form by taking the nonsingular matrix X as 98: * 99: * X = Q*( I 0 ) 100: * ( 0 inv(R) ) 101: * 102: * Arguments 103: * ========= 104: * 105: * JOBU (input) CHARACTER*1 106: * = 'U': Unitary matrix U is computed; 107: * = 'N': U is not computed. 108: * 109: * JOBV (input) CHARACTER*1 110: * = 'V': Unitary matrix V is computed; 111: * = 'N': V is not computed. 112: * 113: * JOBQ (input) CHARACTER*1 114: * = 'Q': Unitary matrix Q is computed; 115: * = 'N': Q is not computed. 116: * 117: * M (input) INTEGER 118: * The number of rows of the matrix A. M >= 0. 119: * 120: * N (input) INTEGER 121: * The number of columns of the matrices A and B. N >= 0. 122: * 123: * P (input) INTEGER 124: * The number of rows of the matrix B. P >= 0. 125: * 126: * K (output) INTEGER 127: * L (output) INTEGER 128: * On exit, K and L specify the dimension of the subblocks 129: * described in Purpose. 130: * K + L = effective numerical rank of (A',B')'. 131: * 132: * A (input/output) COMPLEX*16 array, dimension (LDA,N) 133: * On entry, the M-by-N matrix A. 134: * On exit, A contains the triangular matrix R, or part of R. 135: * See Purpose for details. 136: * 137: * LDA (input) INTEGER 138: * The leading dimension of the array A. LDA >= max(1,M). 139: * 140: * B (input/output) COMPLEX*16 array, dimension (LDB,N) 141: * On entry, the P-by-N matrix B. 142: * On exit, B contains part of the triangular matrix R if 143: * M-K-L < 0. See Purpose for details. 144: * 145: * LDB (input) INTEGER 146: * The leading dimension of the array B. LDB >= max(1,P). 147: * 148: * ALPHA (output) DOUBLE PRECISION array, dimension (N) 149: * BETA (output) DOUBLE PRECISION array, dimension (N) 150: * On exit, ALPHA and BETA contain the generalized singular 151: * value pairs of A and B; 152: * ALPHA(1:K) = 1, 153: * BETA(1:K) = 0, 154: * and if M-K-L >= 0, 155: * ALPHA(K+1:K+L) = C, 156: * BETA(K+1:K+L) = S, 157: * or if M-K-L < 0, 158: * ALPHA(K+1:M)= C, ALPHA(M+1:K+L)= 0 159: * BETA(K+1:M) = S, BETA(M+1:K+L) = 1 160: * and 161: * ALPHA(K+L+1:N) = 0 162: * BETA(K+L+1:N) = 0 163: * 164: * U (output) COMPLEX*16 array, dimension (LDU,M) 165: * If JOBU = 'U', U contains the M-by-M unitary matrix U. 166: * If JOBU = 'N', U is not referenced. 167: * 168: * LDU (input) INTEGER 169: * The leading dimension of the array U. LDU >= max(1,M) if 170: * JOBU = 'U'; LDU >= 1 otherwise. 171: * 172: * V (output) COMPLEX*16 array, dimension (LDV,P) 173: * If JOBV = 'V', V contains the P-by-P unitary matrix V. 174: * If JOBV = 'N', V is not referenced. 175: * 176: * LDV (input) INTEGER 177: * The leading dimension of the array V. LDV >= max(1,P) if 178: * JOBV = 'V'; LDV >= 1 otherwise. 179: * 180: * Q (output) COMPLEX*16 array, dimension (LDQ,N) 181: * If JOBQ = 'Q', Q contains the N-by-N unitary matrix Q. 182: * If JOBQ = 'N', Q is not referenced. 183: * 184: * LDQ (input) INTEGER 185: * The leading dimension of the array Q. LDQ >= max(1,N) if 186: * JOBQ = 'Q'; LDQ >= 1 otherwise. 187: * 188: * WORK (workspace) COMPLEX*16 array, dimension (max(3*N,M,P)+N) 189: * 190: * RWORK (workspace) DOUBLE PRECISION array, dimension (2*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 ZTGSJA. 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, ZLANGE 235: EXTERNAL LSAME, DLAMCH, ZLANGE 236: * .. 237: * .. External Subroutines .. 238: EXTERNAL DCOPY, XERBLA, ZGGSVP, ZTGSJA 239: * .. 240: * .. Intrinsic Functions .. 241: INTRINSIC MAX, MIN 242: * .. 243: * .. Executable Statements .. 244: * 245: * Decode and 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( 'ZGGSVD', -INFO ) 277: RETURN 278: END IF 279: * 280: * Compute the Frobenius norm of matrices A and B 281: * 282: ANORM = ZLANGE( '1', M, N, A, LDA, RWORK ) 283: BNORM = ZLANGE( '1', P, N, B, LDB, RWORK ) 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: CALL ZGGSVP( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, TOLA, 294: $ TOLB, K, L, U, LDU, V, LDV, Q, LDQ, IWORK, RWORK, 295: $ WORK, WORK( N+1 ), INFO ) 296: * 297: * Compute the GSVD of two upper "triangular" matrices 298: * 299: CALL ZTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, LDB, 300: $ TOLA, TOLB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, 301: $ WORK, NCYCLE, INFO ) 302: * 303: * Sort the singular values and store the pivot indices in IWORK 304: * Copy ALPHA to RWORK, then sort ALPHA in RWORK 305: * 306: CALL DCOPY( N, ALPHA, 1, RWORK, 1 ) 307: IBND = MIN( L, M-K ) 308: DO 20 I = 1, IBND 309: * 310: * Scan for largest ALPHA(K+I) 311: * 312: ISUB = I 313: SMAX = RWORK( K+I ) 314: DO 10 J = I + 1, IBND 315: TEMP = RWORK( K+J ) 316: IF( TEMP.GT.SMAX ) THEN 317: ISUB = J 318: SMAX = TEMP 319: END IF 320: 10 CONTINUE 321: IF( ISUB.NE.I ) THEN 322: RWORK( K+ISUB ) = RWORK( K+I ) 323: RWORK( K+I ) = SMAX 324: IWORK( K+I ) = K + ISUB 325: ELSE 326: IWORK( K+I ) = K + I 327: END IF 328: 20 CONTINUE 329: * 330: RETURN 331: * 332: * End of ZGGSVD 333: * 334: END