 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CGGSVD( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, 00002 $ LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, 00003 $ RWORK, IWORK, INFO ) 00004 * 00005 * -- LAPACK driver routine (version 3.3.1) -- 00006 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00007 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00008 * -- April 2011 -- 00009 * 00010 * .. Scalar Arguments .. 00011 CHARACTER JOBQ, JOBU, JOBV 00012 INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P 00013 * .. 00014 * .. Array Arguments .. 00015 INTEGER IWORK( * ) 00016 REAL ALPHA( * ), BETA( * ), RWORK( * ) 00017 COMPLEX A( LDA, * ), B( LDB, * ), Q( LDQ, * ), 00018 $ U( LDU, * ), V( LDV, * ), WORK( * ) 00019 * .. 00020 * 00021 * Purpose 00022 * ======= 00023 * 00024 * CGGSVD computes the generalized singular value decomposition (GSVD) 00025 * of an M-by-N complex matrix A and P-by-N complex matrix B: 00026 * 00027 * U**H*A*Q = D1*( 0 R ), V**H*B*Q = D2*( 0 R ) 00028 * 00029 * where U, V and Q are unitary matrices. 00030 * Let K+L = the effective numerical rank of the 00031 * matrix (A**H,B**H)**H, then R is a (K+L)-by-(K+L) nonsingular upper 00032 * triangular matrix, D1 and D2 are M-by-(K+L) and P-by-(K+L) "diagonal" 00033 * matrices and of the following structures, respectively: 00034 * 00035 * If M-K-L >= 0, 00036 * 00037 * K L 00038 * D1 = K ( I 0 ) 00039 * L ( 0 C ) 00040 * M-K-L ( 0 0 ) 00041 * 00042 * K L 00043 * D2 = L ( 0 S ) 00044 * P-L ( 0 0 ) 00045 * 00046 * N-K-L K L 00047 * ( 0 R ) = K ( 0 R11 R12 ) 00048 * L ( 0 0 R22 ) 00049 * 00050 * where 00051 * 00052 * C = diag( ALPHA(K+1), ... , ALPHA(K+L) ), 00053 * S = diag( BETA(K+1), ... , BETA(K+L) ), 00054 * C**2 + S**2 = I. 00055 * 00056 * R is stored in A(1:K+L,N-K-L+1:N) on exit. 00057 * 00058 * If M-K-L < 0, 00059 * 00060 * K M-K K+L-M 00061 * D1 = K ( I 0 0 ) 00062 * M-K ( 0 C 0 ) 00063 * 00064 * K M-K K+L-M 00065 * D2 = M-K ( 0 S 0 ) 00066 * K+L-M ( 0 0 I ) 00067 * P-L ( 0 0 0 ) 00068 * 00069 * N-K-L K M-K K+L-M 00070 * ( 0 R ) = K ( 0 R11 R12 R13 ) 00071 * M-K ( 0 0 R22 R23 ) 00072 * K+L-M ( 0 0 0 R33 ) 00073 * 00074 * where 00075 * 00076 * C = diag( ALPHA(K+1), ... , ALPHA(M) ), 00077 * S = diag( BETA(K+1), ... , BETA(M) ), 00078 * C**2 + S**2 = I. 00079 * 00080 * (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), and R33 is stored 00081 * ( 0 R22 R23 ) 00082 * in B(M-K+1:L,N+M-K-L+1:N) on exit. 00083 * 00084 * The routine computes C, S, R, and optionally the unitary 00085 * transformation matrices U, V and Q. 00086 * 00087 * In particular, if B is an N-by-N nonsingular matrix, then the GSVD of 00088 * A and B implicitly gives the SVD of A*inv(B): 00089 * A*inv(B) = U*(D1*inv(D2))*V**H. 00090 * If ( A**H,B**H)**H has orthnormal columns, then the GSVD of A and B is also 00091 * equal to the CS decomposition of A and B. Furthermore, the GSVD can 00092 * be used to derive the solution of the eigenvalue problem: 00093 * A**H*A x = lambda* B**H*B x. 00094 * In some literature, the GSVD of A and B is presented in the form 00095 * U**H*A*X = ( 0 D1 ), V**H*B*X = ( 0 D2 ) 00096 * where U and V are orthogonal and X is nonsingular, and D1 and D2 are 00097 * ``diagonal''. The former GSVD form can be converted to the latter 00098 * form by taking the nonsingular matrix X as 00099 * 00100 * X = Q*( I 0 ) 00101 * ( 0 inv(R) ) 00102 * 00103 * Arguments 00104 * ========= 00105 * 00106 * JOBU (input) CHARACTER*1 00107 * = 'U': Unitary matrix U is computed; 00108 * = 'N': U is not computed. 00109 * 00110 * JOBV (input) CHARACTER*1 00111 * = 'V': Unitary matrix V is computed; 00112 * = 'N': V is not computed. 00113 * 00114 * JOBQ (input) CHARACTER*1 00115 * = 'Q': Unitary matrix Q is computed; 00116 * = 'N': Q is not computed. 00117 * 00118 * M (input) INTEGER 00119 * The number of rows of the matrix A. M >= 0. 00120 * 00121 * N (input) INTEGER 00122 * The number of columns of the matrices A and B. N >= 0. 00123 * 00124 * P (input) INTEGER 00125 * The number of rows of the matrix B. P >= 0. 00126 * 00127 * K (output) INTEGER 00128 * L (output) INTEGER 00129 * On exit, K and L specify the dimension of the subblocks 00130 * described in Purpose. 00131 * K + L = effective numerical rank of (A**H,B**H)**H. 00132 * 00133 * A (input/output) COMPLEX array, dimension (LDA,N) 00134 * On entry, the M-by-N matrix A. 00135 * On exit, A contains the triangular matrix R, or part of R. 00136 * See Purpose for details. 00137 * 00138 * LDA (input) INTEGER 00139 * The leading dimension of the array A. LDA >= max(1,M). 00140 * 00141 * B (input/output) COMPLEX array, dimension (LDB,N) 00142 * On entry, the P-by-N matrix B. 00143 * On exit, B contains part of the triangular matrix R if 00144 * M-K-L < 0. See Purpose for details. 00145 * 00146 * LDB (input) INTEGER 00147 * The leading dimension of the array B. LDB >= max(1,P). 00148 * 00149 * ALPHA (output) REAL array, dimension (N) 00150 * BETA (output) REAL array, dimension (N) 00151 * On exit, ALPHA and BETA contain the generalized singular 00152 * value pairs of A and B; 00153 * ALPHA(1:K) = 1, 00154 * BETA(1:K) = 0, 00155 * and if M-K-L >= 0, 00156 * ALPHA(K+1:K+L) = C, 00157 * BETA(K+1:K+L) = S, 00158 * or if M-K-L < 0, 00159 * ALPHA(K+1:M)= C, ALPHA(M+1:K+L)= 0 00160 * BETA(K+1:M) = S, BETA(M+1:K+L) = 1 00161 * and 00162 * ALPHA(K+L+1:N) = 0 00163 * BETA(K+L+1:N) = 0 00164 * 00165 * U (output) COMPLEX array, dimension (LDU,M) 00166 * If JOBU = 'U', U contains the M-by-M unitary matrix U. 00167 * If JOBU = 'N', U is not referenced. 00168 * 00169 * LDU (input) INTEGER 00170 * The leading dimension of the array U. LDU >= max(1,M) if 00171 * JOBU = 'U'; LDU >= 1 otherwise. 00172 * 00173 * V (output) COMPLEX array, dimension (LDV,P) 00174 * If JOBV = 'V', V contains the P-by-P unitary matrix V. 00175 * If JOBV = 'N', V is not referenced. 00176 * 00177 * LDV (input) INTEGER 00178 * The leading dimension of the array V. LDV >= max(1,P) if 00179 * JOBV = 'V'; LDV >= 1 otherwise. 00180 * 00181 * Q (output) COMPLEX array, dimension (LDQ,N) 00182 * If JOBQ = 'Q', Q contains the N-by-N unitary matrix Q. 00183 * If JOBQ = 'N', Q is not referenced. 00184 * 00185 * LDQ (input) INTEGER 00186 * The leading dimension of the array Q. LDQ >= max(1,N) if 00187 * JOBQ = 'Q'; LDQ >= 1 otherwise. 00188 * 00189 * WORK (workspace) COMPLEX array, dimension (max(3*N,M,P)+N) 00190 * 00191 * RWORK (workspace) REAL array, dimension (2*N) 00192 * 00193 * IWORK (workspace/output) INTEGER array, dimension (N) 00194 * On exit, IWORK stores the sorting information. More 00195 * precisely, the following loop will sort ALPHA 00196 * for I = K+1, min(M,K+L) 00197 * swap ALPHA(I) and ALPHA(IWORK(I)) 00198 * endfor 00199 * such that ALPHA(1) >= ALPHA(2) >= ... >= ALPHA(N). 00200 * 00201 * INFO (output) INTEGER 00202 * = 0: successful exit. 00203 * < 0: if INFO = -i, the i-th argument had an illegal value. 00204 * > 0: if INFO = 1, the Jacobi-type procedure failed to 00205 * converge. For further details, see subroutine CTGSJA. 00206 * 00207 * Internal Parameters 00208 * =================== 00209 * 00210 * TOLA REAL 00211 * TOLB REAL 00212 * TOLA and TOLB are the thresholds to determine the effective 00213 * rank of (A**H,B**H)**H. Generally, they are set to 00214 * TOLA = MAX(M,N)*norm(A)*MACHEPS, 00215 * TOLB = MAX(P,N)*norm(B)*MACHEPS. 00216 * The size of TOLA and TOLB may affect the size of backward 00217 * errors of the decomposition. 00218 * 00219 * Further Details 00220 * =============== 00221 * 00222 * 2-96 Based on modifications by 00223 * Ming Gu and Huan Ren, Computer Science Division, University of 00224 * California at Berkeley, USA 00225 * 00226 * ===================================================================== 00227 * 00228 * .. Local Scalars .. 00229 LOGICAL WANTQ, WANTU, WANTV 00230 INTEGER I, IBND, ISUB, J, NCYCLE 00231 REAL ANORM, BNORM, SMAX, TEMP, TOLA, TOLB, ULP, UNFL 00232 * .. 00233 * .. External Functions .. 00234 LOGICAL LSAME 00235 REAL CLANGE, SLAMCH 00236 EXTERNAL LSAME, CLANGE, SLAMCH 00237 * .. 00238 * .. External Subroutines .. 00239 EXTERNAL CGGSVP, CTGSJA, SCOPY, XERBLA 00240 * .. 00241 * .. Intrinsic Functions .. 00242 INTRINSIC MAX, MIN 00243 * .. 00244 * .. Executable Statements .. 00245 * 00246 * Decode and test the input parameters 00247 * 00248 WANTU = LSAME( JOBU, 'U' ) 00249 WANTV = LSAME( JOBV, 'V' ) 00250 WANTQ = LSAME( JOBQ, 'Q' ) 00251 * 00252 INFO = 0 00253 IF( .NOT.( WANTU .OR. LSAME( JOBU, 'N' ) ) ) THEN 00254 INFO = -1 00255 ELSE IF( .NOT.( WANTV .OR. LSAME( JOBV, 'N' ) ) ) THEN 00256 INFO = -2 00257 ELSE IF( .NOT.( WANTQ .OR. LSAME( JOBQ, 'N' ) ) ) THEN 00258 INFO = -3 00259 ELSE IF( M.LT.0 ) THEN 00260 INFO = -4 00261 ELSE IF( N.LT.0 ) THEN 00262 INFO = -5 00263 ELSE IF( P.LT.0 ) THEN 00264 INFO = -6 00265 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00266 INFO = -10 00267 ELSE IF( LDB.LT.MAX( 1, P ) ) THEN 00268 INFO = -12 00269 ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN 00270 INFO = -16 00271 ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN 00272 INFO = -18 00273 ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN 00274 INFO = -20 00275 END IF 00276 IF( INFO.NE.0 ) THEN 00277 CALL XERBLA( 'CGGSVD', -INFO ) 00278 RETURN 00279 END IF 00280 * 00281 * Compute the Frobenius norm of matrices A and B 00282 * 00283 ANORM = CLANGE( '1', M, N, A, LDA, RWORK ) 00284 BNORM = CLANGE( '1', P, N, B, LDB, RWORK ) 00285 * 00286 * Get machine precision and set up threshold for determining 00287 * the effective numerical rank of the matrices A and B. 00288 * 00289 ULP = SLAMCH( 'Precision' ) 00290 UNFL = SLAMCH( 'Safe Minimum' ) 00291 TOLA = MAX( M, N )*MAX( ANORM, UNFL )*ULP 00292 TOLB = MAX( P, N )*MAX( BNORM, UNFL )*ULP 00293 * 00294 CALL CGGSVP( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, TOLA, 00295 $ TOLB, K, L, U, LDU, V, LDV, Q, LDQ, IWORK, RWORK, 00296 $ WORK, WORK( N+1 ), INFO ) 00297 * 00298 * Compute the GSVD of two upper "triangular" matrices 00299 * 00300 CALL CTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, LDB, 00301 $ TOLA, TOLB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, 00302 $ WORK, NCYCLE, INFO ) 00303 * 00304 * Sort the singular values and store the pivot indices in IWORK 00305 * Copy ALPHA to RWORK, then sort ALPHA in RWORK 00306 * 00307 CALL SCOPY( N, ALPHA, 1, RWORK, 1 ) 00308 IBND = MIN( L, M-K ) 00309 DO 20 I = 1, IBND 00310 * 00311 * Scan for largest ALPHA(K+I) 00312 * 00313 ISUB = I 00314 SMAX = RWORK( K+I ) 00315 DO 10 J = I + 1, IBND 00316 TEMP = RWORK( K+J ) 00317 IF( TEMP.GT.SMAX ) THEN 00318 ISUB = J 00319 SMAX = TEMP 00320 END IF 00321 10 CONTINUE 00322 IF( ISUB.NE.I ) THEN 00323 RWORK( K+ISUB ) = RWORK( K+I ) 00324 RWORK( K+I ) = SMAX 00325 IWORK( K+I ) = K + ISUB 00326 ELSE 00327 IWORK( K+I ) = K + I 00328 END IF 00329 20 CONTINUE 00330 * 00331 RETURN 00332 * 00333 * End of CGGSVD 00334 * 00335 END
1.7.3 
