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