![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief <b> CGGSVD computes the singular value decomposition (SVD) for OTHER matrices</b> 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download CGGSVD + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cggsvd.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cggsvd.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cggsvd.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE CGGSVD( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, 00022 * LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, 00023 * RWORK, IWORK, INFO ) 00024 * 00025 * .. Scalar Arguments .. 00026 * CHARACTER JOBQ, JOBU, JOBV 00027 * INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P 00028 * .. 00029 * .. Array Arguments .. 00030 * INTEGER IWORK( * ) 00031 * REAL ALPHA( * ), BETA( * ), RWORK( * ) 00032 * COMPLEX A( LDA, * ), B( LDB, * ), Q( LDQ, * ), 00033 * $ U( LDU, * ), V( LDV, * ), WORK( * ) 00034 * .. 00035 * 00036 * 00037 *> \par Purpose: 00038 * ============= 00039 *> 00040 *> \verbatim 00041 *> 00042 *> CGGSVD computes the generalized singular value decomposition (GSVD) 00043 *> of an M-by-N complex matrix A and P-by-N complex matrix B: 00044 *> 00045 *> U**H*A*Q = D1*( 0 R ), V**H*B*Q = D2*( 0 R ) 00046 *> 00047 *> where U, V and Q are unitary matrices. 00048 *> Let K+L = the effective numerical rank of the 00049 *> matrix (A**H,B**H)**H, then R is a (K+L)-by-(K+L) nonsingular upper 00050 *> triangular matrix, D1 and D2 are M-by-(K+L) and P-by-(K+L) "diagonal" 00051 *> matrices and of the following structures, respectively: 00052 *> 00053 *> If M-K-L >= 0, 00054 *> 00055 *> K L 00056 *> D1 = K ( I 0 ) 00057 *> L ( 0 C ) 00058 *> M-K-L ( 0 0 ) 00059 *> 00060 *> K L 00061 *> D2 = L ( 0 S ) 00062 *> P-L ( 0 0 ) 00063 *> 00064 *> N-K-L K L 00065 *> ( 0 R ) = K ( 0 R11 R12 ) 00066 *> L ( 0 0 R22 ) 00067 *> 00068 *> where 00069 *> 00070 *> C = diag( ALPHA(K+1), ... , ALPHA(K+L) ), 00071 *> S = diag( BETA(K+1), ... , BETA(K+L) ), 00072 *> C**2 + S**2 = I. 00073 *> 00074 *> R is stored in A(1:K+L,N-K-L+1:N) on exit. 00075 *> 00076 *> If M-K-L < 0, 00077 *> 00078 *> K M-K K+L-M 00079 *> D1 = K ( I 0 0 ) 00080 *> M-K ( 0 C 0 ) 00081 *> 00082 *> K M-K K+L-M 00083 *> D2 = M-K ( 0 S 0 ) 00084 *> K+L-M ( 0 0 I ) 00085 *> P-L ( 0 0 0 ) 00086 *> 00087 *> N-K-L K M-K K+L-M 00088 *> ( 0 R ) = K ( 0 R11 R12 R13 ) 00089 *> M-K ( 0 0 R22 R23 ) 00090 *> K+L-M ( 0 0 0 R33 ) 00091 *> 00092 *> where 00093 *> 00094 *> C = diag( ALPHA(K+1), ... , ALPHA(M) ), 00095 *> S = diag( BETA(K+1), ... , BETA(M) ), 00096 *> C**2 + S**2 = I. 00097 *> 00098 *> (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), and R33 is stored 00099 *> ( 0 R22 R23 ) 00100 *> in B(M-K+1:L,N+M-K-L+1:N) on exit. 00101 *> 00102 *> The routine computes C, S, R, and optionally the unitary 00103 *> transformation matrices U, V and Q. 00104 *> 00105 *> In particular, if B is an N-by-N nonsingular matrix, then the GSVD of 00106 *> A and B implicitly gives the SVD of A*inv(B): 00107 *> A*inv(B) = U*(D1*inv(D2))*V**H. 00108 *> If ( A**H,B**H)**H has orthnormal columns, then the GSVD of A and B is also 00109 *> equal to the CS decomposition of A and B. Furthermore, the GSVD can 00110 *> be used to derive the solution of the eigenvalue problem: 00111 *> A**H*A x = lambda* B**H*B x. 00112 *> In some literature, the GSVD of A and B is presented in the form 00113 *> U**H*A*X = ( 0 D1 ), V**H*B*X = ( 0 D2 ) 00114 *> where U and V are orthogonal and X is nonsingular, and D1 and D2 are 00115 *> ``diagonal''. The former GSVD form can be converted to the latter 00116 *> form by taking the nonsingular matrix X as 00117 *> 00118 *> X = Q*( I 0 ) 00119 *> ( 0 inv(R) ) 00120 *> \endverbatim 00121 * 00122 * Arguments: 00123 * ========== 00124 * 00125 *> \param[in] JOBU 00126 *> \verbatim 00127 *> JOBU is CHARACTER*1 00128 *> = 'U': Unitary matrix U is computed; 00129 *> = 'N': U is not computed. 00130 *> \endverbatim 00131 *> 00132 *> \param[in] JOBV 00133 *> \verbatim 00134 *> JOBV is CHARACTER*1 00135 *> = 'V': Unitary matrix V is computed; 00136 *> = 'N': V is not computed. 00137 *> \endverbatim 00138 *> 00139 *> \param[in] JOBQ 00140 *> \verbatim 00141 *> JOBQ is CHARACTER*1 00142 *> = 'Q': Unitary matrix Q is computed; 00143 *> = 'N': Q is not computed. 00144 *> \endverbatim 00145 *> 00146 *> \param[in] M 00147 *> \verbatim 00148 *> M is INTEGER 00149 *> The number of rows of the matrix A. M >= 0. 00150 *> \endverbatim 00151 *> 00152 *> \param[in] N 00153 *> \verbatim 00154 *> N is INTEGER 00155 *> The number of columns of the matrices A and B. N >= 0. 00156 *> \endverbatim 00157 *> 00158 *> \param[in] P 00159 *> \verbatim 00160 *> P is INTEGER 00161 *> The number of rows of the matrix B. P >= 0. 00162 *> \endverbatim 00163 *> 00164 *> \param[out] K 00165 *> \verbatim 00166 *> K is INTEGER 00167 *> \endverbatim 00168 *> 00169 *> \param[out] L 00170 *> \verbatim 00171 *> L is INTEGER 00172 *> 00173 *> On exit, K and L specify the dimension of the subblocks 00174 *> described in Purpose. 00175 *> K + L = effective numerical rank of (A**H,B**H)**H. 00176 *> \endverbatim 00177 *> 00178 *> \param[in,out] A 00179 *> \verbatim 00180 *> A is COMPLEX array, dimension (LDA,N) 00181 *> On entry, the M-by-N matrix A. 00182 *> On exit, A contains the triangular matrix R, or part of R. 00183 *> See Purpose for details. 00184 *> \endverbatim 00185 *> 00186 *> \param[in] LDA 00187 *> \verbatim 00188 *> LDA is INTEGER 00189 *> The leading dimension of the array A. LDA >= max(1,M). 00190 *> \endverbatim 00191 *> 00192 *> \param[in,out] B 00193 *> \verbatim 00194 *> B is COMPLEX array, dimension (LDB,N) 00195 *> On entry, the P-by-N matrix B. 00196 *> On exit, B contains part of the triangular matrix R if 00197 *> M-K-L < 0. See Purpose for details. 00198 *> \endverbatim 00199 *> 00200 *> \param[in] LDB 00201 *> \verbatim 00202 *> LDB is INTEGER 00203 *> The leading dimension of the array B. LDB >= max(1,P). 00204 *> \endverbatim 00205 *> 00206 *> \param[out] ALPHA 00207 *> \verbatim 00208 *> ALPHA is REAL array, dimension (N) 00209 *> \endverbatim 00210 *> 00211 *> \param[out] BETA 00212 *> \verbatim 00213 *> BETA is REAL array, dimension (N) 00214 *> 00215 *> On exit, ALPHA and BETA contain the generalized singular 00216 *> value pairs of A and B; 00217 *> ALPHA(1:K) = 1, 00218 *> BETA(1:K) = 0, 00219 *> and if M-K-L >= 0, 00220 *> ALPHA(K+1:K+L) = C, 00221 *> BETA(K+1:K+L) = S, 00222 *> or if M-K-L < 0, 00223 *> ALPHA(K+1:M)=C, ALPHA(M+1:K+L)=0 00224 *> BETA(K+1:M) =S, BETA(M+1:K+L) =1 00225 *> and 00226 *> ALPHA(K+L+1:N) = 0 00227 *> BETA(K+L+1:N) = 0 00228 *> \endverbatim 00229 *> 00230 *> \param[out] U 00231 *> \verbatim 00232 *> U is COMPLEX array, dimension (LDU,M) 00233 *> If JOBU = 'U', U contains the M-by-M unitary matrix U. 00234 *> If JOBU = 'N', U is not referenced. 00235 *> \endverbatim 00236 *> 00237 *> \param[in] LDU 00238 *> \verbatim 00239 *> LDU is INTEGER 00240 *> The leading dimension of the array U. LDU >= max(1,M) if 00241 *> JOBU = 'U'; LDU >= 1 otherwise. 00242 *> \endverbatim 00243 *> 00244 *> \param[out] V 00245 *> \verbatim 00246 *> V is COMPLEX array, dimension (LDV,P) 00247 *> If JOBV = 'V', V contains the P-by-P unitary matrix V. 00248 *> If JOBV = 'N', V is not referenced. 00249 *> \endverbatim 00250 *> 00251 *> \param[in] LDV 00252 *> \verbatim 00253 *> LDV is INTEGER 00254 *> The leading dimension of the array V. LDV >= max(1,P) if 00255 *> JOBV = 'V'; LDV >= 1 otherwise. 00256 *> \endverbatim 00257 *> 00258 *> \param[out] Q 00259 *> \verbatim 00260 *> Q is COMPLEX array, dimension (LDQ,N) 00261 *> If JOBQ = 'Q', Q contains the N-by-N unitary matrix Q. 00262 *> If JOBQ = 'N', Q is not referenced. 00263 *> \endverbatim 00264 *> 00265 *> \param[in] LDQ 00266 *> \verbatim 00267 *> LDQ is INTEGER 00268 *> The leading dimension of the array Q. LDQ >= max(1,N) if 00269 *> JOBQ = 'Q'; LDQ >= 1 otherwise. 00270 *> \endverbatim 00271 *> 00272 *> \param[out] WORK 00273 *> \verbatim 00274 *> WORK is COMPLEX array, dimension (max(3*N,M,P)+N) 00275 *> \endverbatim 00276 *> 00277 *> \param[out] RWORK 00278 *> \verbatim 00279 *> RWORK is REAL array, dimension (2*N) 00280 *> \endverbatim 00281 *> 00282 *> \param[out] IWORK 00283 *> \verbatim 00284 *> IWORK is INTEGER array, dimension (N) 00285 *> On exit, IWORK stores the sorting information. More 00286 *> precisely, the following loop will sort ALPHA 00287 *> for I = K+1, min(M,K+L) 00288 *> swap ALPHA(I) and ALPHA(IWORK(I)) 00289 *> endfor 00290 *> such that ALPHA(1) >= ALPHA(2) >= ... >= ALPHA(N). 00291 *> \endverbatim 00292 *> 00293 *> \param[out] INFO 00294 *> \verbatim 00295 *> INFO is INTEGER 00296 *> = 0: successful exit. 00297 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00298 *> > 0: if INFO = 1, the Jacobi-type procedure failed to 00299 *> converge. For further details, see subroutine CTGSJA. 00300 *> \endverbatim 00301 * 00302 *> \par Internal Parameters: 00303 * ========================= 00304 *> 00305 *> \verbatim 00306 *> TOLA REAL 00307 *> TOLB REAL 00308 *> TOLA and TOLB are the thresholds to determine the effective 00309 *> rank of (A**H,B**H)**H. Generally, they are set to 00310 *> TOLA = MAX(M,N)*norm(A)*MACHEPS, 00311 *> TOLB = MAX(P,N)*norm(B)*MACHEPS. 00312 *> The size of TOLA and TOLB may affect the size of backward 00313 *> errors of the decomposition. 00314 *> \endverbatim 00315 * 00316 * Authors: 00317 * ======== 00318 * 00319 *> \author Univ. of Tennessee 00320 *> \author Univ. of California Berkeley 00321 *> \author Univ. of Colorado Denver 00322 *> \author NAG Ltd. 00323 * 00324 *> \date November 2011 00325 * 00326 *> \ingroup complexOTHERsing 00327 * 00328 *> \par Contributors: 00329 * ================== 00330 *> 00331 *> Ming Gu and Huan Ren, Computer Science Division, University of 00332 *> California at Berkeley, USA 00333 *> 00334 * ===================================================================== 00335 SUBROUTINE CGGSVD( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, 00336 $ LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, 00337 $ RWORK, IWORK, INFO ) 00338 * 00339 * -- LAPACK driver routine (version 3.4.0) -- 00340 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00341 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00342 * November 2011 00343 * 00344 * .. Scalar Arguments .. 00345 CHARACTER JOBQ, JOBU, JOBV 00346 INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P 00347 * .. 00348 * .. Array Arguments .. 00349 INTEGER IWORK( * ) 00350 REAL ALPHA( * ), BETA( * ), RWORK( * ) 00351 COMPLEX A( LDA, * ), B( LDB, * ), Q( LDQ, * ), 00352 $ U( LDU, * ), V( LDV, * ), WORK( * ) 00353 * .. 00354 * 00355 * ===================================================================== 00356 * 00357 * .. Local Scalars .. 00358 LOGICAL WANTQ, WANTU, WANTV 00359 INTEGER I, IBND, ISUB, J, NCYCLE 00360 REAL ANORM, BNORM, SMAX, TEMP, TOLA, TOLB, ULP, UNFL 00361 * .. 00362 * .. External Functions .. 00363 LOGICAL LSAME 00364 REAL CLANGE, SLAMCH 00365 EXTERNAL LSAME, CLANGE, SLAMCH 00366 * .. 00367 * .. External Subroutines .. 00368 EXTERNAL CGGSVP, CTGSJA, SCOPY, XERBLA 00369 * .. 00370 * .. Intrinsic Functions .. 00371 INTRINSIC MAX, MIN 00372 * .. 00373 * .. Executable Statements .. 00374 * 00375 * Decode and test the input parameters 00376 * 00377 WANTU = LSAME( JOBU, 'U' ) 00378 WANTV = LSAME( JOBV, 'V' ) 00379 WANTQ = LSAME( JOBQ, 'Q' ) 00380 * 00381 INFO = 0 00382 IF( .NOT.( WANTU .OR. LSAME( JOBU, 'N' ) ) ) THEN 00383 INFO = -1 00384 ELSE IF( .NOT.( WANTV .OR. LSAME( JOBV, 'N' ) ) ) THEN 00385 INFO = -2 00386 ELSE IF( .NOT.( WANTQ .OR. LSAME( JOBQ, 'N' ) ) ) THEN 00387 INFO = -3 00388 ELSE IF( M.LT.0 ) THEN 00389 INFO = -4 00390 ELSE IF( N.LT.0 ) THEN 00391 INFO = -5 00392 ELSE IF( P.LT.0 ) THEN 00393 INFO = -6 00394 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00395 INFO = -10 00396 ELSE IF( LDB.LT.MAX( 1, P ) ) THEN 00397 INFO = -12 00398 ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN 00399 INFO = -16 00400 ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN 00401 INFO = -18 00402 ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN 00403 INFO = -20 00404 END IF 00405 IF( INFO.NE.0 ) THEN 00406 CALL XERBLA( 'CGGSVD', -INFO ) 00407 RETURN 00408 END IF 00409 * 00410 * Compute the Frobenius norm of matrices A and B 00411 * 00412 ANORM = CLANGE( '1', M, N, A, LDA, RWORK ) 00413 BNORM = CLANGE( '1', P, N, B, LDB, RWORK ) 00414 * 00415 * Get machine precision and set up threshold for determining 00416 * the effective numerical rank of the matrices A and B. 00417 * 00418 ULP = SLAMCH( 'Precision' ) 00419 UNFL = SLAMCH( 'Safe Minimum' ) 00420 TOLA = MAX( M, N )*MAX( ANORM, UNFL )*ULP 00421 TOLB = MAX( P, N )*MAX( BNORM, UNFL )*ULP 00422 * 00423 CALL CGGSVP( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, TOLA, 00424 $ TOLB, K, L, U, LDU, V, LDV, Q, LDQ, IWORK, RWORK, 00425 $ WORK, WORK( N+1 ), INFO ) 00426 * 00427 * Compute the GSVD of two upper "triangular" matrices 00428 * 00429 CALL CTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, LDB, 00430 $ TOLA, TOLB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, 00431 $ WORK, NCYCLE, INFO ) 00432 * 00433 * Sort the singular values and store the pivot indices in IWORK 00434 * Copy ALPHA to RWORK, then sort ALPHA in RWORK 00435 * 00436 CALL SCOPY( N, ALPHA, 1, RWORK, 1 ) 00437 IBND = MIN( L, M-K ) 00438 DO 20 I = 1, IBND 00439 * 00440 * Scan for largest ALPHA(K+I) 00441 * 00442 ISUB = I 00443 SMAX = RWORK( K+I ) 00444 DO 10 J = I + 1, IBND 00445 TEMP = RWORK( K+J ) 00446 IF( TEMP.GT.SMAX ) THEN 00447 ISUB = J 00448 SMAX = TEMP 00449 END IF 00450 10 CONTINUE 00451 IF( ISUB.NE.I ) THEN 00452 RWORK( K+ISUB ) = RWORK( K+I ) 00453 RWORK( K+I ) = SMAX 00454 IWORK( K+I ) = K + ISUB 00455 ELSE 00456 IWORK( K+I ) = K + I 00457 END IF 00458 20 CONTINUE 00459 * 00460 RETURN 00461 * 00462 * End of CGGSVD 00463 * 00464 END