![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief <b> DGGSVD 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 DGGSVD + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dggsvd.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dggsvd.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dggsvd.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DGGSVD( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, 00022 * LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, 00023 * 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 A( LDA, * ), ALPHA( * ), B( LDB, * ), 00032 * $ BETA( * ), Q( LDQ, * ), U( LDU, * ), 00033 * $ V( LDV, * ), WORK( * ) 00034 * .. 00035 * 00036 * 00037 *> \par Purpose: 00038 * ============= 00039 *> 00040 *> \verbatim 00041 *> 00042 *> DGGSVD computes the generalized singular value decomposition (GSVD) 00043 *> of an M-by-N real matrix A and P-by-N real matrix B: 00044 *> 00045 *> U**T*A*Q = D1*( 0 R ), V**T*B*Q = D2*( 0 R ) 00046 *> 00047 *> where U, V and Q are orthogonal matrices. 00048 *> Let K+L = the effective numerical rank of the matrix (A**T,B**T)**T, 00049 *> then R is a K+L-by-K+L nonsingular upper triangular matrix, D1 and 00050 *> D2 are M-by-(K+L) and P-by-(K+L) "diagonal" matrices and of the 00051 *> 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 orthogonal 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**T. 00108 *> If ( A**T,B**T)**T has orthonormal columns, then the GSVD of A and B is 00109 *> also equal to the CS decomposition of A and B. Furthermore, the GSVD 00110 *> can be used to derive the solution of the eigenvalue problem: 00111 *> A**T*A x = lambda* B**T*B x. 00112 *> In some literature, the GSVD of A and B is presented in the form 00113 *> U**T*A*X = ( 0 D1 ), V**T*B*X = ( 0 D2 ) 00114 *> where U and V are orthogonal and X is nonsingular, 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': Orthogonal 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': Orthogonal 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': Orthogonal 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**T,B**T)**T. 00176 *> \endverbatim 00177 *> 00178 *> \param[in,out] A 00179 *> \verbatim 00180 *> A is DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDB,N) 00195 *> On entry, the P-by-N matrix B. 00196 *> On exit, B contains the triangular matrix R if M-K-L < 0. 00197 *> 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 DOUBLE PRECISION array, dimension (N) 00209 *> \endverbatim 00210 *> 00211 *> \param[out] BETA 00212 *> \verbatim 00213 *> BETA is DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDU,M) 00233 *> If JOBU = 'U', U contains the M-by-M orthogonal 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 DOUBLE PRECISION array, dimension (LDV,P) 00247 *> If JOBV = 'V', V contains the P-by-P orthogonal 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 DOUBLE PRECISION array, dimension (LDQ,N) 00261 *> If JOBQ = 'Q', Q contains the N-by-N orthogonal 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 DOUBLE PRECISION array, 00275 *> dimension (max(3*N,M,P)+N) 00276 *> \endverbatim 00277 *> 00278 *> \param[out] IWORK 00279 *> \verbatim 00280 *> IWORK is INTEGER array, dimension (N) 00281 *> On exit, IWORK stores the sorting information. More 00282 *> precisely, the following loop will sort ALPHA 00283 *> for I = K+1, min(M,K+L) 00284 *> swap ALPHA(I) and ALPHA(IWORK(I)) 00285 *> endfor 00286 *> such that ALPHA(1) >= ALPHA(2) >= ... >= ALPHA(N). 00287 *> \endverbatim 00288 *> 00289 *> \param[out] INFO 00290 *> \verbatim 00291 *> INFO is INTEGER 00292 *> = 0: successful exit 00293 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00294 *> > 0: if INFO = 1, the Jacobi-type procedure failed to 00295 *> converge. For further details, see subroutine DTGSJA. 00296 *> \endverbatim 00297 * 00298 *> \par Internal Parameters: 00299 * ========================= 00300 *> 00301 *> \verbatim 00302 *> TOLA DOUBLE PRECISION 00303 *> TOLB DOUBLE PRECISION 00304 *> TOLA and TOLB are the thresholds to determine the effective 00305 *> rank of (A',B')**T. Generally, they are set to 00306 *> TOLA = MAX(M,N)*norm(A)*MAZHEPS, 00307 *> TOLB = MAX(P,N)*norm(B)*MAZHEPS. 00308 *> The size of TOLA and TOLB may affect the size of backward 00309 *> errors of the decomposition. 00310 *> \endverbatim 00311 * 00312 * Authors: 00313 * ======== 00314 * 00315 *> \author Univ. of Tennessee 00316 *> \author Univ. of California Berkeley 00317 *> \author Univ. of Colorado Denver 00318 *> \author NAG Ltd. 00319 * 00320 *> \date November 2011 00321 * 00322 *> \ingroup doubleOTHERsing 00323 * 00324 *> \par Contributors: 00325 * ================== 00326 *> 00327 *> Ming Gu and Huan Ren, Computer Science Division, University of 00328 *> California at Berkeley, USA 00329 *> 00330 * ===================================================================== 00331 SUBROUTINE DGGSVD( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, 00332 $ LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, 00333 $ IWORK, INFO ) 00334 * 00335 * -- LAPACK driver routine (version 3.4.0) -- 00336 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00337 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00338 * November 2011 00339 * 00340 * .. Scalar Arguments .. 00341 CHARACTER JOBQ, JOBU, JOBV 00342 INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P 00343 * .. 00344 * .. Array Arguments .. 00345 INTEGER IWORK( * ) 00346 DOUBLE PRECISION A( LDA, * ), ALPHA( * ), B( LDB, * ), 00347 $ BETA( * ), Q( LDQ, * ), U( LDU, * ), 00348 $ V( LDV, * ), WORK( * ) 00349 * .. 00350 * 00351 * ===================================================================== 00352 * 00353 * .. Local Scalars .. 00354 LOGICAL WANTQ, WANTU, WANTV 00355 INTEGER I, IBND, ISUB, J, NCYCLE 00356 DOUBLE PRECISION ANORM, BNORM, SMAX, TEMP, TOLA, TOLB, ULP, UNFL 00357 * .. 00358 * .. External Functions .. 00359 LOGICAL LSAME 00360 DOUBLE PRECISION DLAMCH, DLANGE 00361 EXTERNAL LSAME, DLAMCH, DLANGE 00362 * .. 00363 * .. External Subroutines .. 00364 EXTERNAL DCOPY, DGGSVP, DTGSJA, XERBLA 00365 * .. 00366 * .. Intrinsic Functions .. 00367 INTRINSIC MAX, MIN 00368 * .. 00369 * .. Executable Statements .. 00370 * 00371 * Test the input parameters 00372 * 00373 WANTU = LSAME( JOBU, 'U' ) 00374 WANTV = LSAME( JOBV, 'V' ) 00375 WANTQ = LSAME( JOBQ, 'Q' ) 00376 * 00377 INFO = 0 00378 IF( .NOT.( WANTU .OR. LSAME( JOBU, 'N' ) ) ) THEN 00379 INFO = -1 00380 ELSE IF( .NOT.( WANTV .OR. LSAME( JOBV, 'N' ) ) ) THEN 00381 INFO = -2 00382 ELSE IF( .NOT.( WANTQ .OR. LSAME( JOBQ, 'N' ) ) ) THEN 00383 INFO = -3 00384 ELSE IF( M.LT.0 ) THEN 00385 INFO = -4 00386 ELSE IF( N.LT.0 ) THEN 00387 INFO = -5 00388 ELSE IF( P.LT.0 ) THEN 00389 INFO = -6 00390 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00391 INFO = -10 00392 ELSE IF( LDB.LT.MAX( 1, P ) ) THEN 00393 INFO = -12 00394 ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN 00395 INFO = -16 00396 ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN 00397 INFO = -18 00398 ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN 00399 INFO = -20 00400 END IF 00401 IF( INFO.NE.0 ) THEN 00402 CALL XERBLA( 'DGGSVD', -INFO ) 00403 RETURN 00404 END IF 00405 * 00406 * Compute the Frobenius norm of matrices A and B 00407 * 00408 ANORM = DLANGE( '1', M, N, A, LDA, WORK ) 00409 BNORM = DLANGE( '1', P, N, B, LDB, WORK ) 00410 * 00411 * Get machine precision and set up threshold for determining 00412 * the effective numerical rank of the matrices A and B. 00413 * 00414 ULP = DLAMCH( 'Precision' ) 00415 UNFL = DLAMCH( 'Safe Minimum' ) 00416 TOLA = MAX( M, N )*MAX( ANORM, UNFL )*ULP 00417 TOLB = MAX( P, N )*MAX( BNORM, UNFL )*ULP 00418 * 00419 * Preprocessing 00420 * 00421 CALL DGGSVP( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, TOLA, 00422 $ TOLB, K, L, U, LDU, V, LDV, Q, LDQ, IWORK, WORK, 00423 $ WORK( N+1 ), INFO ) 00424 * 00425 * Compute the GSVD of two upper "triangular" matrices 00426 * 00427 CALL DTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, LDB, 00428 $ TOLA, TOLB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, 00429 $ WORK, NCYCLE, INFO ) 00430 * 00431 * Sort the singular values and store the pivot indices in IWORK 00432 * Copy ALPHA to WORK, then sort ALPHA in WORK 00433 * 00434 CALL DCOPY( N, ALPHA, 1, WORK, 1 ) 00435 IBND = MIN( L, M-K ) 00436 DO 20 I = 1, IBND 00437 * 00438 * Scan for largest ALPHA(K+I) 00439 * 00440 ISUB = I 00441 SMAX = WORK( K+I ) 00442 DO 10 J = I + 1, IBND 00443 TEMP = WORK( K+J ) 00444 IF( TEMP.GT.SMAX ) THEN 00445 ISUB = J 00446 SMAX = TEMP 00447 END IF 00448 10 CONTINUE 00449 IF( ISUB.NE.I ) THEN 00450 WORK( K+ISUB ) = WORK( K+I ) 00451 WORK( K+I ) = SMAX 00452 IWORK( K+I ) = K + ISUB 00453 ELSE 00454 IWORK( K+I ) = K + I 00455 END IF 00456 20 CONTINUE 00457 * 00458 RETURN 00459 * 00460 * End of DGGSVD 00461 * 00462 END