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