00001 SUBROUTINE CTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, 00002 $ LDB, TOLA, TOLB, ALPHA, BETA, U, LDU, V, LDV, 00003 $ Q, LDQ, WORK, NCYCLE, INFO ) 00004 * 00005 * -- LAPACK routine (version 3.2.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 2009 -- 00009 * 00010 * .. Scalar Arguments .. 00011 CHARACTER JOBQ, JOBU, JOBV 00012 INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, 00013 $ NCYCLE, P 00014 REAL TOLA, TOLB 00015 * .. 00016 * .. Array Arguments .. 00017 REAL ALPHA( * ), BETA( * ) 00018 COMPLEX A( LDA, * ), B( LDB, * ), Q( LDQ, * ), 00019 $ U( LDU, * ), V( LDV, * ), WORK( * ) 00020 * .. 00021 * 00022 * Purpose 00023 * ======= 00024 * 00025 * CTGSJA computes the generalized singular value decomposition (GSVD) 00026 * of two complex upper triangular (or trapezoidal) matrices A and B. 00027 * 00028 * On entry, it is assumed that matrices A and B have the following 00029 * forms, which may be obtained by the preprocessing subroutine CGGSVP 00030 * from a general M-by-N matrix A and P-by-N matrix B: 00031 * 00032 * N-K-L K L 00033 * A = K ( 0 A12 A13 ) if M-K-L >= 0; 00034 * L ( 0 0 A23 ) 00035 * M-K-L ( 0 0 0 ) 00036 * 00037 * N-K-L K L 00038 * A = K ( 0 A12 A13 ) if M-K-L < 0; 00039 * M-K ( 0 0 A23 ) 00040 * 00041 * N-K-L K L 00042 * B = L ( 0 0 B13 ) 00043 * P-L ( 0 0 0 ) 00044 * 00045 * where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular 00046 * upper triangular; A23 is L-by-L upper triangular if M-K-L >= 0, 00047 * otherwise A23 is (M-K)-by-L upper trapezoidal. 00048 * 00049 * On exit, 00050 * 00051 * U'*A*Q = D1*( 0 R ), V'*B*Q = D2*( 0 R ), 00052 * 00053 * where U, V and Q are unitary matrices, Z' denotes the conjugate 00054 * transpose of Z, R is a nonsingular upper triangular matrix, and D1 00055 * and D2 are ``diagonal'' matrices, which are of the following 00056 * structures: 00057 * 00058 * If M-K-L >= 0, 00059 * 00060 * K L 00061 * D1 = K ( I 0 ) 00062 * L ( 0 C ) 00063 * M-K-L ( 0 0 ) 00064 * 00065 * K L 00066 * D2 = L ( 0 S ) 00067 * P-L ( 0 0 ) 00068 * 00069 * N-K-L K L 00070 * ( 0 R ) = K ( 0 R11 R12 ) K 00071 * L ( 0 0 R22 ) L 00072 * 00073 * where 00074 * 00075 * C = diag( ALPHA(K+1), ... , ALPHA(K+L) ), 00076 * S = diag( BETA(K+1), ... , BETA(K+L) ), 00077 * C**2 + S**2 = I. 00078 * 00079 * R is stored in A(1:K+L,N-K-L+1:N) on exit. 00080 * 00081 * If M-K-L < 0, 00082 * 00083 * K M-K K+L-M 00084 * D1 = K ( I 0 0 ) 00085 * M-K ( 0 C 0 ) 00086 * 00087 * K M-K K+L-M 00088 * D2 = M-K ( 0 S 0 ) 00089 * K+L-M ( 0 0 I ) 00090 * P-L ( 0 0 0 ) 00091 * 00092 * N-K-L K M-K K+L-M 00093 * ( 0 R ) = K ( 0 R11 R12 R13 ) 00094 * M-K ( 0 0 R22 R23 ) 00095 * K+L-M ( 0 0 0 R33 ) 00096 * 00097 * where 00098 * C = diag( ALPHA(K+1), ... , ALPHA(M) ), 00099 * S = diag( BETA(K+1), ... , BETA(M) ), 00100 * C**2 + S**2 = I. 00101 * 00102 * R = ( R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N) and R33 is stored 00103 * ( 0 R22 R23 ) 00104 * in B(M-K+1:L,N+M-K-L+1:N) on exit. 00105 * 00106 * The computation of the unitary transformation matrices U, V or Q 00107 * is optional. These matrices may either be formed explicitly, or they 00108 * may be postmultiplied into input matrices U1, V1, or Q1. 00109 * 00110 * Arguments 00111 * ========= 00112 * 00113 * JOBU (input) CHARACTER*1 00114 * = 'U': U must contain a unitary matrix U1 on entry, and 00115 * the product U1*U is returned; 00116 * = 'I': U is initialized to the unit matrix, and the 00117 * unitary matrix U is returned; 00118 * = 'N': U is not computed. 00119 * 00120 * JOBV (input) CHARACTER*1 00121 * = 'V': V must contain a unitary matrix V1 on entry, and 00122 * the product V1*V is returned; 00123 * = 'I': V is initialized to the unit matrix, and the 00124 * unitary matrix V is returned; 00125 * = 'N': V is not computed. 00126 * 00127 * JOBQ (input) CHARACTER*1 00128 * = 'Q': Q must contain a unitary matrix Q1 on entry, and 00129 * the product Q1*Q is returned; 00130 * = 'I': Q is initialized to the unit matrix, and the 00131 * unitary matrix Q is returned; 00132 * = 'N': Q is not computed. 00133 * 00134 * M (input) INTEGER 00135 * The number of rows of the matrix A. M >= 0. 00136 * 00137 * P (input) INTEGER 00138 * The number of rows of the matrix B. P >= 0. 00139 * 00140 * N (input) INTEGER 00141 * The number of columns of the matrices A and B. N >= 0. 00142 * 00143 * K (input) INTEGER 00144 * L (input) INTEGER 00145 * K and L specify the subblocks in the input matrices A and B: 00146 * A23 = A(K+1:MIN(K+L,M),N-L+1:N) and B13 = B(1:L,,N-L+1:N) 00147 * of A and B, whose GSVD is going to be computed by CTGSJA. 00148 * See Further Details. 00149 * 00150 * A (input/output) COMPLEX array, dimension (LDA,N) 00151 * On entry, the M-by-N matrix A. 00152 * On exit, A(N-K+1:N,1:MIN(K+L,M) ) contains the triangular 00153 * matrix R or part of R. See Purpose for details. 00154 * 00155 * LDA (input) INTEGER 00156 * The leading dimension of the array A. LDA >= max(1,M). 00157 * 00158 * B (input/output) COMPLEX array, dimension (LDB,N) 00159 * On entry, the P-by-N matrix B. 00160 * On exit, if necessary, B(M-K+1:L,N+M-K-L+1:N) contains 00161 * a part of R. See Purpose for details. 00162 * 00163 * LDB (input) INTEGER 00164 * The leading dimension of the array B. LDB >= max(1,P). 00165 * 00166 * TOLA (input) REAL 00167 * TOLB (input) REAL 00168 * TOLA and TOLB are the convergence criteria for the Jacobi- 00169 * Kogbetliantz iteration procedure. Generally, they are the 00170 * same as used in the preprocessing step, say 00171 * TOLA = MAX(M,N)*norm(A)*MACHEPS, 00172 * TOLB = MAX(P,N)*norm(B)*MACHEPS. 00173 * 00174 * ALPHA (output) REAL array, dimension (N) 00175 * BETA (output) REAL array, dimension (N) 00176 * On exit, ALPHA and BETA contain the generalized singular 00177 * value pairs of A and B; 00178 * ALPHA(1:K) = 1, 00179 * BETA(1:K) = 0, 00180 * and if M-K-L >= 0, 00181 * ALPHA(K+1:K+L) = diag(C), 00182 * BETA(K+1:K+L) = diag(S), 00183 * or if M-K-L < 0, 00184 * ALPHA(K+1:M)= C, ALPHA(M+1:K+L)= 0 00185 * BETA(K+1:M) = S, BETA(M+1:K+L) = 1. 00186 * Furthermore, if K+L < N, 00187 * ALPHA(K+L+1:N) = 0 00188 * BETA(K+L+1:N) = 0. 00189 * 00190 * U (input/output) COMPLEX array, dimension (LDU,M) 00191 * On entry, if JOBU = 'U', U must contain a matrix U1 (usually 00192 * the unitary matrix returned by CGGSVP). 00193 * On exit, 00194 * if JOBU = 'I', U contains the unitary matrix U; 00195 * if JOBU = 'U', U contains the product U1*U. 00196 * If JOBU = 'N', U is not referenced. 00197 * 00198 * LDU (input) INTEGER 00199 * The leading dimension of the array U. LDU >= max(1,M) if 00200 * JOBU = 'U'; LDU >= 1 otherwise. 00201 * 00202 * V (input/output) COMPLEX array, dimension (LDV,P) 00203 * On entry, if JOBV = 'V', V must contain a matrix V1 (usually 00204 * the unitary matrix returned by CGGSVP). 00205 * On exit, 00206 * if JOBV = 'I', V contains the unitary matrix V; 00207 * if JOBV = 'V', V contains the product V1*V. 00208 * If JOBV = 'N', V is not referenced. 00209 * 00210 * LDV (input) INTEGER 00211 * The leading dimension of the array V. LDV >= max(1,P) if 00212 * JOBV = 'V'; LDV >= 1 otherwise. 00213 * 00214 * Q (input/output) COMPLEX array, dimension (LDQ,N) 00215 * On entry, if JOBQ = 'Q', Q must contain a matrix Q1 (usually 00216 * the unitary matrix returned by CGGSVP). 00217 * On exit, 00218 * if JOBQ = 'I', Q contains the unitary matrix Q; 00219 * if JOBQ = 'Q', Q contains the product Q1*Q. 00220 * If JOBQ = 'N', Q is not referenced. 00221 * 00222 * LDQ (input) INTEGER 00223 * The leading dimension of the array Q. LDQ >= max(1,N) if 00224 * JOBQ = 'Q'; LDQ >= 1 otherwise. 00225 * 00226 * WORK (workspace) COMPLEX array, dimension (2*N) 00227 * 00228 * NCYCLE (output) INTEGER 00229 * The number of cycles required for convergence. 00230 * 00231 * INFO (output) INTEGER 00232 * = 0: successful exit 00233 * < 0: if INFO = -i, the i-th argument had an illegal value. 00234 * = 1: the procedure does not converge after MAXIT cycles. 00235 * 00236 * Internal Parameters 00237 * =================== 00238 * 00239 * MAXIT INTEGER 00240 * MAXIT specifies the total loops that the iterative procedure 00241 * may take. If after MAXIT cycles, the routine fails to 00242 * converge, we return INFO = 1. 00243 * 00244 * Further Details 00245 * =============== 00246 * 00247 * CTGSJA essentially uses a variant of Kogbetliantz algorithm to reduce 00248 * min(L,M-K)-by-L triangular (or trapezoidal) matrix A23 and L-by-L 00249 * matrix B13 to the form: 00250 * 00251 * U1'*A13*Q1 = C1*R1; V1'*B13*Q1 = S1*R1, 00252 * 00253 * where U1, V1 and Q1 are unitary matrix, and Z' is the conjugate 00254 * transpose of Z. C1 and S1 are diagonal matrices satisfying 00255 * 00256 * C1**2 + S1**2 = I, 00257 * 00258 * and R1 is an L-by-L nonsingular upper triangular matrix. 00259 * 00260 * ===================================================================== 00261 * 00262 * .. Parameters .. 00263 INTEGER MAXIT 00264 PARAMETER ( MAXIT = 40 ) 00265 REAL ZERO, ONE 00266 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00267 COMPLEX CZERO, CONE 00268 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 00269 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 00270 * .. 00271 * .. Local Scalars .. 00272 * 00273 LOGICAL INITQ, INITU, INITV, UPPER, WANTQ, WANTU, WANTV 00274 INTEGER I, J, KCYCLE 00275 REAL A1, A3, B1, B3, CSQ, CSU, CSV, ERROR, GAMMA, 00276 $ RWK, SSMIN 00277 COMPLEX A2, B2, SNQ, SNU, SNV 00278 * .. 00279 * .. External Functions .. 00280 LOGICAL LSAME 00281 EXTERNAL LSAME 00282 * .. 00283 * .. External Subroutines .. 00284 EXTERNAL CCOPY, CLAGS2, CLAPLL, CLASET, CROT, CSSCAL, 00285 $ SLARTG, XERBLA 00286 * .. 00287 * .. Intrinsic Functions .. 00288 INTRINSIC ABS, CONJG, MAX, MIN, REAL 00289 * .. 00290 * .. Executable Statements .. 00291 * 00292 * Decode and test the input parameters 00293 * 00294 INITU = LSAME( JOBU, 'I' ) 00295 WANTU = INITU .OR. LSAME( JOBU, 'U' ) 00296 * 00297 INITV = LSAME( JOBV, 'I' ) 00298 WANTV = INITV .OR. LSAME( JOBV, 'V' ) 00299 * 00300 INITQ = LSAME( JOBQ, 'I' ) 00301 WANTQ = INITQ .OR. LSAME( JOBQ, 'Q' ) 00302 * 00303 INFO = 0 00304 IF( .NOT.( INITU .OR. WANTU .OR. LSAME( JOBU, 'N' ) ) ) THEN 00305 INFO = -1 00306 ELSE IF( .NOT.( INITV .OR. WANTV .OR. LSAME( JOBV, 'N' ) ) ) THEN 00307 INFO = -2 00308 ELSE IF( .NOT.( INITQ .OR. WANTQ .OR. LSAME( JOBQ, 'N' ) ) ) THEN 00309 INFO = -3 00310 ELSE IF( M.LT.0 ) THEN 00311 INFO = -4 00312 ELSE IF( P.LT.0 ) THEN 00313 INFO = -5 00314 ELSE IF( N.LT.0 ) THEN 00315 INFO = -6 00316 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00317 INFO = -10 00318 ELSE IF( LDB.LT.MAX( 1, P ) ) THEN 00319 INFO = -12 00320 ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN 00321 INFO = -18 00322 ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN 00323 INFO = -20 00324 ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN 00325 INFO = -22 00326 END IF 00327 IF( INFO.NE.0 ) THEN 00328 CALL XERBLA( 'CTGSJA', -INFO ) 00329 RETURN 00330 END IF 00331 * 00332 * Initialize U, V and Q, if necessary 00333 * 00334 IF( INITU ) 00335 $ CALL CLASET( 'Full', M, M, CZERO, CONE, U, LDU ) 00336 IF( INITV ) 00337 $ CALL CLASET( 'Full', P, P, CZERO, CONE, V, LDV ) 00338 IF( INITQ ) 00339 $ CALL CLASET( 'Full', N, N, CZERO, CONE, Q, LDQ ) 00340 * 00341 * Loop until convergence 00342 * 00343 UPPER = .FALSE. 00344 DO 40 KCYCLE = 1, MAXIT 00345 * 00346 UPPER = .NOT.UPPER 00347 * 00348 DO 20 I = 1, L - 1 00349 DO 10 J = I + 1, L 00350 * 00351 A1 = ZERO 00352 A2 = CZERO 00353 A3 = ZERO 00354 IF( K+I.LE.M ) 00355 $ A1 = REAL( A( K+I, N-L+I ) ) 00356 IF( K+J.LE.M ) 00357 $ A3 = REAL( A( K+J, N-L+J ) ) 00358 * 00359 B1 = REAL( B( I, N-L+I ) ) 00360 B3 = REAL( B( J, N-L+J ) ) 00361 * 00362 IF( UPPER ) THEN 00363 IF( K+I.LE.M ) 00364 $ A2 = A( K+I, N-L+J ) 00365 B2 = B( I, N-L+J ) 00366 ELSE 00367 IF( K+J.LE.M ) 00368 $ A2 = A( K+J, N-L+I ) 00369 B2 = B( J, N-L+I ) 00370 END IF 00371 * 00372 CALL CLAGS2( UPPER, A1, A2, A3, B1, B2, B3, CSU, SNU, 00373 $ CSV, SNV, CSQ, SNQ ) 00374 * 00375 * Update (K+I)-th and (K+J)-th rows of matrix A: U'*A 00376 * 00377 IF( K+J.LE.M ) 00378 $ CALL CROT( L, A( K+J, N-L+1 ), LDA, A( K+I, N-L+1 ), 00379 $ LDA, CSU, CONJG( SNU ) ) 00380 * 00381 * Update I-th and J-th rows of matrix B: V'*B 00382 * 00383 CALL CROT( L, B( J, N-L+1 ), LDB, B( I, N-L+1 ), LDB, 00384 $ CSV, CONJG( SNV ) ) 00385 * 00386 * Update (N-L+I)-th and (N-L+J)-th columns of matrices 00387 * A and B: A*Q and B*Q 00388 * 00389 CALL CROT( MIN( K+L, M ), A( 1, N-L+J ), 1, 00390 $ A( 1, N-L+I ), 1, CSQ, SNQ ) 00391 * 00392 CALL CROT( L, B( 1, N-L+J ), 1, B( 1, N-L+I ), 1, CSQ, 00393 $ SNQ ) 00394 * 00395 IF( UPPER ) THEN 00396 IF( K+I.LE.M ) 00397 $ A( K+I, N-L+J ) = CZERO 00398 B( I, N-L+J ) = CZERO 00399 ELSE 00400 IF( K+J.LE.M ) 00401 $ A( K+J, N-L+I ) = CZERO 00402 B( J, N-L+I ) = CZERO 00403 END IF 00404 * 00405 * Ensure that the diagonal elements of A and B are real. 00406 * 00407 IF( K+I.LE.M ) 00408 $ A( K+I, N-L+I ) = REAL( A( K+I, N-L+I ) ) 00409 IF( K+J.LE.M ) 00410 $ A( K+J, N-L+J ) = REAL( A( K+J, N-L+J ) ) 00411 B( I, N-L+I ) = REAL( B( I, N-L+I ) ) 00412 B( J, N-L+J ) = REAL( B( J, N-L+J ) ) 00413 * 00414 * Update unitary matrices U, V, Q, if desired. 00415 * 00416 IF( WANTU .AND. K+J.LE.M ) 00417 $ CALL CROT( M, U( 1, K+J ), 1, U( 1, K+I ), 1, CSU, 00418 $ SNU ) 00419 * 00420 IF( WANTV ) 00421 $ CALL CROT( P, V( 1, J ), 1, V( 1, I ), 1, CSV, SNV ) 00422 * 00423 IF( WANTQ ) 00424 $ CALL CROT( N, Q( 1, N-L+J ), 1, Q( 1, N-L+I ), 1, CSQ, 00425 $ SNQ ) 00426 * 00427 10 CONTINUE 00428 20 CONTINUE 00429 * 00430 IF( .NOT.UPPER ) THEN 00431 * 00432 * The matrices A13 and B13 were lower triangular at the start 00433 * of the cycle, and are now upper triangular. 00434 * 00435 * Convergence test: test the parallelism of the corresponding 00436 * rows of A and B. 00437 * 00438 ERROR = ZERO 00439 DO 30 I = 1, MIN( L, M-K ) 00440 CALL CCOPY( L-I+1, A( K+I, N-L+I ), LDA, WORK, 1 ) 00441 CALL CCOPY( L-I+1, B( I, N-L+I ), LDB, WORK( L+1 ), 1 ) 00442 CALL CLAPLL( L-I+1, WORK, 1, WORK( L+1 ), 1, SSMIN ) 00443 ERROR = MAX( ERROR, SSMIN ) 00444 30 CONTINUE 00445 * 00446 IF( ABS( ERROR ).LE.MIN( TOLA, TOLB ) ) 00447 $ GO TO 50 00448 END IF 00449 * 00450 * End of cycle loop 00451 * 00452 40 CONTINUE 00453 * 00454 * The algorithm has not converged after MAXIT cycles. 00455 * 00456 INFO = 1 00457 GO TO 100 00458 * 00459 50 CONTINUE 00460 * 00461 * If ERROR <= MIN(TOLA,TOLB), then the algorithm has converged. 00462 * Compute the generalized singular value pairs (ALPHA, BETA), and 00463 * set the triangular matrix R to array A. 00464 * 00465 DO 60 I = 1, K 00466 ALPHA( I ) = ONE 00467 BETA( I ) = ZERO 00468 60 CONTINUE 00469 * 00470 DO 70 I = 1, MIN( L, M-K ) 00471 * 00472 A1 = REAL( A( K+I, N-L+I ) ) 00473 B1 = REAL( B( I, N-L+I ) ) 00474 * 00475 IF( A1.NE.ZERO ) THEN 00476 GAMMA = B1 / A1 00477 * 00478 IF( GAMMA.LT.ZERO ) THEN 00479 CALL CSSCAL( L-I+1, -ONE, B( I, N-L+I ), LDB ) 00480 IF( WANTV ) 00481 $ CALL CSSCAL( P, -ONE, V( 1, I ), 1 ) 00482 END IF 00483 * 00484 CALL SLARTG( ABS( GAMMA ), ONE, BETA( K+I ), ALPHA( K+I ), 00485 $ RWK ) 00486 * 00487 IF( ALPHA( K+I ).GE.BETA( K+I ) ) THEN 00488 CALL CSSCAL( L-I+1, ONE / ALPHA( K+I ), A( K+I, N-L+I ), 00489 $ LDA ) 00490 ELSE 00491 CALL CSSCAL( L-I+1, ONE / BETA( K+I ), B( I, N-L+I ), 00492 $ LDB ) 00493 CALL CCOPY( L-I+1, B( I, N-L+I ), LDB, A( K+I, N-L+I ), 00494 $ LDA ) 00495 END IF 00496 * 00497 ELSE 00498 ALPHA( K+I ) = ZERO 00499 BETA( K+I ) = ONE 00500 CALL CCOPY( L-I+1, B( I, N-L+I ), LDB, A( K+I, N-L+I ), 00501 $ LDA ) 00502 END IF 00503 70 CONTINUE 00504 * 00505 * Post-assignment 00506 * 00507 DO 80 I = M + 1, K + L 00508 ALPHA( I ) = ZERO 00509 BETA( I ) = ONE 00510 80 CONTINUE 00511 * 00512 IF( K+L.LT.N ) THEN 00513 DO 90 I = K + L + 1, N 00514 ALPHA( I ) = ZERO 00515 BETA( I ) = ZERO 00516 90 CONTINUE 00517 END IF 00518 * 00519 100 CONTINUE 00520 NCYCLE = KCYCLE 00521 * 00522 RETURN 00523 * 00524 * End of CTGSJA 00525 * 00526 END
1.7.2