| 1 | SUBROUTINE CGGSVD( 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.1) --
|
|---|
| 6 | * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
|
|---|
| 7 | * November 2006
|
|---|
| 8 | *
|
|---|
| 9 | * .. Scalar Arguments ..
|
|---|
| 10 | CHARACTER JOBQ, JOBU, JOBV
|
|---|
| 11 | INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P
|
|---|
| 12 | * ..
|
|---|
| 13 | * .. Array Arguments ..
|
|---|
| 14 | INTEGER IWORK( * )
|
|---|
| 15 | REAL ALPHA( * ), BETA( * ), RWORK( * )
|
|---|
| 16 | COMPLEX A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
|
|---|
| 17 | $ U( LDU, * ), V( LDV, * ), WORK( * )
|
|---|
| 18 | * ..
|
|---|
| 19 | *
|
|---|
| 20 | * Purpose
|
|---|
| 21 | * =======
|
|---|
| 22 | *
|
|---|
| 23 | * CGGSVD computes the generalized singular value decomposition (GSVD)
|
|---|
| 24 | * of an M-by-N complex matrix A and P-by-N complex matrix B:
|
|---|
| 25 | *
|
|---|
| 26 | * U'*A*Q = D1*( 0 R ), V'*B*Q = D2*( 0 R )
|
|---|
| 27 | *
|
|---|
| 28 | * where U, V and Q are unitary matrices, and Z' means the conjugate
|
|---|
| 29 | * transpose of Z. Let K+L = the effective numerical rank of the
|
|---|
| 30 | * matrix (A',B')', then R is a (K+L)-by-(K+L) nonsingular upper
|
|---|
| 31 | * triangular matrix, D1 and D2 are M-by-(K+L) and P-by-(K+L) "diagonal"
|
|---|
| 32 | * matrices and of the following structures, respectively:
|
|---|
| 33 | *
|
|---|
| 34 | * If M-K-L >= 0,
|
|---|
| 35 | *
|
|---|
| 36 | * K L
|
|---|
| 37 | * D1 = K ( I 0 )
|
|---|
| 38 | * L ( 0 C )
|
|---|
| 39 | * M-K-L ( 0 0 )
|
|---|
| 40 | *
|
|---|
| 41 | * K L
|
|---|
| 42 | * D2 = L ( 0 S )
|
|---|
| 43 | * P-L ( 0 0 )
|
|---|
| 44 | *
|
|---|
| 45 | * N-K-L K L
|
|---|
| 46 | * ( 0 R ) = K ( 0 R11 R12 )
|
|---|
| 47 | * L ( 0 0 R22 )
|
|---|
| 48 | * where
|
|---|
| 49 | *
|
|---|
| 50 | * C = diag( ALPHA(K+1), ... , ALPHA(K+L) ),
|
|---|
| 51 | * S = diag( BETA(K+1), ... , BETA(K+L) ),
|
|---|
| 52 | * C**2 + S**2 = I.
|
|---|
| 53 | *
|
|---|
| 54 | * R is stored in A(1:K+L,N-K-L+1:N) on exit.
|
|---|
| 55 | *
|
|---|
| 56 | * If M-K-L < 0,
|
|---|
| 57 | *
|
|---|
| 58 | * K M-K K+L-M
|
|---|
| 59 | * D1 = K ( I 0 0 )
|
|---|
| 60 | * M-K ( 0 C 0 )
|
|---|
| 61 | *
|
|---|
| 62 | * K M-K K+L-M
|
|---|
| 63 | * D2 = M-K ( 0 S 0 )
|
|---|
| 64 | * K+L-M ( 0 0 I )
|
|---|
| 65 | * P-L ( 0 0 0 )
|
|---|
| 66 | *
|
|---|
| 67 | * N-K-L K M-K K+L-M
|
|---|
| 68 | * ( 0 R ) = K ( 0 R11 R12 R13 )
|
|---|
| 69 | * M-K ( 0 0 R22 R23 )
|
|---|
| 70 | * K+L-M ( 0 0 0 R33 )
|
|---|
| 71 | *
|
|---|
| 72 | * where
|
|---|
| 73 | *
|
|---|
| 74 | * C = diag( ALPHA(K+1), ... , ALPHA(M) ),
|
|---|
| 75 | * S = diag( BETA(K+1), ... , BETA(M) ),
|
|---|
| 76 | * C**2 + S**2 = I.
|
|---|
| 77 | *
|
|---|
| 78 | * (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), and R33 is stored
|
|---|
| 79 | * ( 0 R22 R23 )
|
|---|
| 80 | * in B(M-K+1:L,N+M-K-L+1:N) on exit.
|
|---|
| 81 | *
|
|---|
| 82 | * The routine computes C, S, R, and optionally the unitary
|
|---|
| 83 | * transformation matrices U, V and Q.
|
|---|
| 84 | *
|
|---|
| 85 | * In particular, if B is an N-by-N nonsingular matrix, then the GSVD of
|
|---|
| 86 | * A and B implicitly gives the SVD of A*inv(B):
|
|---|
| 87 | * A*inv(B) = U*(D1*inv(D2))*V'.
|
|---|
| 88 | * If ( A',B')' has orthnormal columns, then the GSVD of A and B is also
|
|---|
| 89 | * equal to the CS decomposition of A and B. Furthermore, the GSVD can
|
|---|
| 90 | * be used to derive the solution of the eigenvalue problem:
|
|---|
| 91 | * A'*A x = lambda* B'*B x.
|
|---|
| 92 | * In some literature, the GSVD of A and B is presented in the form
|
|---|
| 93 | * U'*A*X = ( 0 D1 ), V'*B*X = ( 0 D2 )
|
|---|
| 94 | * where U and V are orthogonal and X is nonsingular, and D1 and D2 are
|
|---|
| 95 | * ``diagonal''. The former GSVD form can be converted to the latter
|
|---|
| 96 | * form by taking the nonsingular matrix X as
|
|---|
| 97 | *
|
|---|
| 98 | * X = Q*( I 0 )
|
|---|
| 99 | * ( 0 inv(R) )
|
|---|
| 100 | *
|
|---|
| 101 | * Arguments
|
|---|
| 102 | * =========
|
|---|
| 103 | *
|
|---|
| 104 | * JOBU (input) CHARACTER*1
|
|---|
| 105 | * = 'U': Unitary matrix U is computed;
|
|---|
| 106 | * = 'N': U is not computed.
|
|---|
| 107 | *
|
|---|
| 108 | * JOBV (input) CHARACTER*1
|
|---|
| 109 | * = 'V': Unitary matrix V is computed;
|
|---|
| 110 | * = 'N': V is not computed.
|
|---|
| 111 | *
|
|---|
| 112 | * JOBQ (input) CHARACTER*1
|
|---|
| 113 | * = 'Q': Unitary matrix Q is computed;
|
|---|
| 114 | * = 'N': Q is not computed.
|
|---|
| 115 | *
|
|---|
| 116 | * M (input) INTEGER
|
|---|
| 117 | * The number of rows of the matrix A. M >= 0.
|
|---|
| 118 | *
|
|---|
| 119 | * N (input) INTEGER
|
|---|
| 120 | * The number of columns of the matrices A and B. N >= 0.
|
|---|
| 121 | *
|
|---|
| 122 | * P (input) INTEGER
|
|---|
| 123 | * The number of rows of the matrix B. P >= 0.
|
|---|
| 124 | *
|
|---|
| 125 | * K (output) INTEGER
|
|---|
| 126 | * L (output) INTEGER
|
|---|
| 127 | * On exit, K and L specify the dimension of the subblocks
|
|---|
| 128 | * described in Purpose.
|
|---|
| 129 | * K + L = effective numerical rank of (A',B')'.
|
|---|
| 130 | *
|
|---|
| 131 | * A (input/output) COMPLEX array, dimension (LDA,N)
|
|---|
| 132 | * On entry, the M-by-N matrix A.
|
|---|
| 133 | * On exit, A contains the triangular matrix R, or part of R.
|
|---|
| 134 | * See Purpose for details.
|
|---|
| 135 | *
|
|---|
| 136 | * LDA (input) INTEGER
|
|---|
| 137 | * The leading dimension of the array A. LDA >= max(1,M).
|
|---|
| 138 | *
|
|---|
| 139 | * B (input/output) COMPLEX array, dimension (LDB,N)
|
|---|
| 140 | * On entry, the P-by-N matrix B.
|
|---|
| 141 | * On exit, B contains part of the triangular matrix R if
|
|---|
| 142 | * M-K-L < 0. See Purpose for details.
|
|---|
| 143 | *
|
|---|
| 144 | * LDB (input) INTEGER
|
|---|
| 145 | * The leading dimension of the array B. LDB >= max(1,P).
|
|---|
| 146 | *
|
|---|
| 147 | * ALPHA (output) REAL array, dimension (N)
|
|---|
| 148 | * BETA (output) REAL array, dimension (N)
|
|---|
| 149 | * On exit, ALPHA and BETA contain the generalized singular
|
|---|
| 150 | * value pairs of A and B;
|
|---|
| 151 | * ALPHA(1:K) = 1,
|
|---|
| 152 | * BETA(1:K) = 0,
|
|---|
| 153 | * and if M-K-L >= 0,
|
|---|
| 154 | * ALPHA(K+1:K+L) = C,
|
|---|
| 155 | * BETA(K+1:K+L) = S,
|
|---|
| 156 | * or if M-K-L < 0,
|
|---|
| 157 | * ALPHA(K+1:M)= C, ALPHA(M+1:K+L)= 0
|
|---|
| 158 | * BETA(K+1:M) = S, BETA(M+1:K+L) = 1
|
|---|
| 159 | * and
|
|---|
| 160 | * ALPHA(K+L+1:N) = 0
|
|---|
| 161 | * BETA(K+L+1:N) = 0
|
|---|
| 162 | *
|
|---|
| 163 | * U (output) COMPLEX array, dimension (LDU,M)
|
|---|
| 164 | * If JOBU = 'U', U contains the M-by-M unitary matrix U.
|
|---|
| 165 | * If JOBU = 'N', U is not referenced.
|
|---|
| 166 | *
|
|---|
| 167 | * LDU (input) INTEGER
|
|---|
| 168 | * The leading dimension of the array U. LDU >= max(1,M) if
|
|---|
| 169 | * JOBU = 'U'; LDU >= 1 otherwise.
|
|---|
| 170 | *
|
|---|
| 171 | * V (output) COMPLEX array, dimension (LDV,P)
|
|---|
| 172 | * If JOBV = 'V', V contains the P-by-P unitary matrix V.
|
|---|
| 173 | * If JOBV = 'N', V is not referenced.
|
|---|
| 174 | *
|
|---|
| 175 | * LDV (input) INTEGER
|
|---|
| 176 | * The leading dimension of the array V. LDV >= max(1,P) if
|
|---|
| 177 | * JOBV = 'V'; LDV >= 1 otherwise.
|
|---|
| 178 | *
|
|---|
| 179 | * Q (output) COMPLEX array, dimension (LDQ,N)
|
|---|
| 180 | * If JOBQ = 'Q', Q contains the N-by-N unitary matrix Q.
|
|---|
| 181 | * If JOBQ = 'N', Q is not referenced.
|
|---|
| 182 | *
|
|---|
| 183 | * LDQ (input) INTEGER
|
|---|
| 184 | * The leading dimension of the array Q. LDQ >= max(1,N) if
|
|---|
| 185 | * JOBQ = 'Q'; LDQ >= 1 otherwise.
|
|---|
| 186 | *
|
|---|
| 187 | * WORK (workspace) COMPLEX array, dimension (max(3*N,M,P)+N)
|
|---|
| 188 | *
|
|---|
| 189 | * RWORK (workspace) REAL array, dimension (2*N)
|
|---|
| 190 | *
|
|---|
| 191 | * IWORK (workspace/output) INTEGER array, dimension (N)
|
|---|
| 192 | * On exit, IWORK stores the sorting information. More
|
|---|
| 193 | * precisely, the following loop will sort ALPHA
|
|---|
| 194 | * for I = K+1, min(M,K+L)
|
|---|
| 195 | * swap ALPHA(I) and ALPHA(IWORK(I))
|
|---|
| 196 | * endfor
|
|---|
| 197 | * such that ALPHA(1) >= ALPHA(2) >= ... >= ALPHA(N).
|
|---|
| 198 | *
|
|---|
| 199 | * INFO (output) INTEGER
|
|---|
| 200 | * = 0: successful exit.
|
|---|
| 201 | * < 0: if INFO = -i, the i-th argument had an illegal value.
|
|---|
| 202 | * > 0: if INFO = 1, the Jacobi-type procedure failed to
|
|---|
| 203 | * converge. For further details, see subroutine CTGSJA.
|
|---|
| 204 | *
|
|---|
| 205 | * Internal Parameters
|
|---|
| 206 | * ===================
|
|---|
| 207 | *
|
|---|
| 208 | * TOLA REAL
|
|---|
| 209 | * TOLB REAL
|
|---|
| 210 | * TOLA and TOLB are the thresholds to determine the effective
|
|---|
| 211 | * rank of (A',B')'. Generally, they are set to
|
|---|
| 212 | * TOLA = MAX(M,N)*norm(A)*MACHEPS,
|
|---|
| 213 | * TOLB = MAX(P,N)*norm(B)*MACHEPS.
|
|---|
| 214 | * The size of TOLA and TOLB may affect the size of backward
|
|---|
| 215 | * errors of the decomposition.
|
|---|
| 216 | *
|
|---|
| 217 | * Further Details
|
|---|
| 218 | * ===============
|
|---|
| 219 | *
|
|---|
| 220 | * 2-96 Based on modifications by
|
|---|
| 221 | * Ming Gu and Huan Ren, Computer Science Division, University of
|
|---|
| 222 | * California at Berkeley, USA
|
|---|
| 223 | *
|
|---|
| 224 | * =====================================================================
|
|---|
| 225 | *
|
|---|
| 226 | * .. Local Scalars ..
|
|---|
| 227 | LOGICAL WANTQ, WANTU, WANTV
|
|---|
| 228 | INTEGER I, IBND, ISUB, J, NCYCLE
|
|---|
| 229 | REAL ANORM, BNORM, SMAX, TEMP, TOLA, TOLB, ULP, UNFL
|
|---|
| 230 | * ..
|
|---|
| 231 | * .. External Functions ..
|
|---|
| 232 | LOGICAL LSAME
|
|---|
| 233 | REAL CLANGE, SLAMCH
|
|---|
| 234 | EXTERNAL LSAME, CLANGE, SLAMCH
|
|---|
| 235 | * ..
|
|---|
| 236 | * .. External Subroutines ..
|
|---|
| 237 | EXTERNAL CGGSVP, CTGSJA, SCOPY, XERBLA
|
|---|
| 238 | * ..
|
|---|
| 239 | * .. Intrinsic Functions ..
|
|---|
| 240 | INTRINSIC MAX, MIN
|
|---|
| 241 | * ..
|
|---|
| 242 | * .. Executable Statements ..
|
|---|
| 243 | *
|
|---|
| 244 | * Decode and test the input parameters
|
|---|
| 245 | *
|
|---|
| 246 | WANTU = LSAME( JOBU, 'U' )
|
|---|
| 247 | WANTV = LSAME( JOBV, 'V' )
|
|---|
| 248 | WANTQ = LSAME( JOBQ, 'Q' )
|
|---|
| 249 | *
|
|---|
| 250 | INFO = 0
|
|---|
| 251 | IF( .NOT.( WANTU .OR. LSAME( JOBU, 'N' ) ) ) THEN
|
|---|
| 252 | INFO = -1
|
|---|
| 253 | ELSE IF( .NOT.( WANTV .OR. LSAME( JOBV, 'N' ) ) ) THEN
|
|---|
| 254 | INFO = -2
|
|---|
| 255 | ELSE IF( .NOT.( WANTQ .OR. LSAME( JOBQ, 'N' ) ) ) THEN
|
|---|
| 256 | INFO = -3
|
|---|
| 257 | ELSE IF( M.LT.0 ) THEN
|
|---|
| 258 | INFO = -4
|
|---|
| 259 | ELSE IF( N.LT.0 ) THEN
|
|---|
| 260 | INFO = -5
|
|---|
| 261 | ELSE IF( P.LT.0 ) THEN
|
|---|
| 262 | INFO = -6
|
|---|
| 263 | ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
|
|---|
| 264 | INFO = -10
|
|---|
| 265 | ELSE IF( LDB.LT.MAX( 1, P ) ) THEN
|
|---|
| 266 | INFO = -12
|
|---|
| 267 | ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN
|
|---|
| 268 | INFO = -16
|
|---|
| 269 | ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN
|
|---|
| 270 | INFO = -18
|
|---|
| 271 | ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
|
|---|
| 272 | INFO = -20
|
|---|
| 273 | END IF
|
|---|
| 274 | IF( INFO.NE.0 ) THEN
|
|---|
| 275 | CALL XERBLA( 'CGGSVD', -INFO )
|
|---|
| 276 | RETURN
|
|---|
| 277 | END IF
|
|---|
| 278 | *
|
|---|
| 279 | * Compute the Frobenius norm of matrices A and B
|
|---|
| 280 | *
|
|---|
| 281 | ANORM = CLANGE( '1', M, N, A, LDA, RWORK )
|
|---|
| 282 | BNORM = CLANGE( '1', P, N, B, LDB, RWORK )
|
|---|
| 283 | *
|
|---|
| 284 | * Get machine precision and set up threshold for determining
|
|---|
| 285 | * the effective numerical rank of the matrices A and B.
|
|---|
| 286 | *
|
|---|
| 287 | ULP = SLAMCH( 'Precision' )
|
|---|
| 288 | UNFL = SLAMCH( 'Safe Minimum' )
|
|---|
| 289 | TOLA = MAX( M, N )*MAX( ANORM, UNFL )*ULP
|
|---|
| 290 | TOLB = MAX( P, N )*MAX( BNORM, UNFL )*ULP
|
|---|
| 291 | *
|
|---|
| 292 | CALL CGGSVP( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, TOLA,
|
|---|
| 293 | $ TOLB, K, L, U, LDU, V, LDV, Q, LDQ, IWORK, RWORK,
|
|---|
| 294 | $ WORK, WORK( N+1 ), INFO )
|
|---|
| 295 | *
|
|---|
| 296 | * Compute the GSVD of two upper "triangular" matrices
|
|---|
| 297 | *
|
|---|
| 298 | CALL CTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, LDB,
|
|---|
| 299 | $ TOLA, TOLB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ,
|
|---|
| 300 | $ WORK, NCYCLE, INFO )
|
|---|
| 301 | *
|
|---|
| 302 | * Sort the singular values and store the pivot indices in IWORK
|
|---|
| 303 | * Copy ALPHA to RWORK, then sort ALPHA in RWORK
|
|---|
| 304 | *
|
|---|
| 305 | CALL SCOPY( N, ALPHA, 1, RWORK, 1 )
|
|---|
| 306 | IBND = MIN( L, M-K )
|
|---|
| 307 | DO 20 I = 1, IBND
|
|---|
| 308 | *
|
|---|
| 309 | * Scan for largest ALPHA(K+I)
|
|---|
| 310 | *
|
|---|
| 311 | ISUB = I
|
|---|
| 312 | SMAX = RWORK( K+I )
|
|---|
| 313 | DO 10 J = I + 1, IBND
|
|---|
| 314 | TEMP = RWORK( K+J )
|
|---|
| 315 | IF( TEMP.GT.SMAX ) THEN
|
|---|
| 316 | ISUB = J
|
|---|
| 317 | SMAX = TEMP
|
|---|
| 318 | END IF
|
|---|
| 319 | 10 CONTINUE
|
|---|
| 320 | IF( ISUB.NE.I ) THEN
|
|---|
| 321 | RWORK( K+ISUB ) = RWORK( K+I )
|
|---|
| 322 | RWORK( K+I ) = SMAX
|
|---|
| 323 | IWORK( K+I ) = K + ISUB
|
|---|
| 324 | ELSE
|
|---|
| 325 | IWORK( K+I ) = K + I
|
|---|
| 326 | END IF
|
|---|
| 327 | 20 CONTINUE
|
|---|
| 328 | *
|
|---|
| 329 | RETURN
|
|---|
| 330 | *
|
|---|
| 331 | * End of CGGSVD
|
|---|
| 332 | *
|
|---|
| 333 | END
|
|---|
![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
