|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE CGGSVP( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, 00002 $ TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ, 00003 $ IWORK, RWORK, TAU, WORK, INFO ) 00004 * 00005 * -- LAPACK routine (version 3.2) -- 00006 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00007 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00008 * November 2006 00009 * 00010 * .. Scalar Arguments .. 00011 CHARACTER JOBQ, JOBU, JOBV 00012 INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P 00013 REAL TOLA, TOLB 00014 * .. 00015 * .. Array Arguments .. 00016 INTEGER IWORK( * ) 00017 REAL RWORK( * ) 00018 COMPLEX A( LDA, * ), B( LDB, * ), Q( LDQ, * ), 00019 $ TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * ) 00020 * .. 00021 * 00022 * Purpose 00023 * ======= 00024 * 00025 * CGGSVP computes unitary matrices U, V and Q such that 00026 * 00027 * N-K-L K L 00028 * U'*A*Q = K ( 0 A12 A13 ) if M-K-L >= 0; 00029 * L ( 0 0 A23 ) 00030 * M-K-L ( 0 0 0 ) 00031 * 00032 * N-K-L K L 00033 * = K ( 0 A12 A13 ) if M-K-L < 0; 00034 * M-K ( 0 0 A23 ) 00035 * 00036 * N-K-L K L 00037 * V'*B*Q = L ( 0 0 B13 ) 00038 * P-L ( 0 0 0 ) 00039 * 00040 * where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular 00041 * upper triangular; A23 is L-by-L upper triangular if M-K-L >= 0, 00042 * otherwise A23 is (M-K)-by-L upper trapezoidal. K+L = the effective 00043 * numerical rank of the (M+P)-by-N matrix (A',B')'. Z' denotes the 00044 * conjugate transpose of Z. 00045 * 00046 * This decomposition is the preprocessing step for computing the 00047 * Generalized Singular Value Decomposition (GSVD), see subroutine 00048 * CGGSVD. 00049 * 00050 * Arguments 00051 * ========= 00052 * 00053 * JOBU (input) CHARACTER*1 00054 * = 'U': Unitary matrix U is computed; 00055 * = 'N': U is not computed. 00056 * 00057 * JOBV (input) CHARACTER*1 00058 * = 'V': Unitary matrix V is computed; 00059 * = 'N': V is not computed. 00060 * 00061 * JOBQ (input) CHARACTER*1 00062 * = 'Q': Unitary matrix Q is computed; 00063 * = 'N': Q is not computed. 00064 * 00065 * M (input) INTEGER 00066 * The number of rows of the matrix A. M >= 0. 00067 * 00068 * P (input) INTEGER 00069 * The number of rows of the matrix B. P >= 0. 00070 * 00071 * N (input) INTEGER 00072 * The number of columns of the matrices A and B. N >= 0. 00073 * 00074 * A (input/output) COMPLEX array, dimension (LDA,N) 00075 * On entry, the M-by-N matrix A. 00076 * On exit, A contains the triangular (or trapezoidal) matrix 00077 * described in the Purpose section. 00078 * 00079 * LDA (input) INTEGER 00080 * The leading dimension of the array A. LDA >= max(1,M). 00081 * 00082 * B (input/output) COMPLEX array, dimension (LDB,N) 00083 * On entry, the P-by-N matrix B. 00084 * On exit, B contains the triangular matrix described in 00085 * the Purpose section. 00086 * 00087 * LDB (input) INTEGER 00088 * The leading dimension of the array B. LDB >= max(1,P). 00089 * 00090 * TOLA (input) REAL 00091 * TOLB (input) REAL 00092 * TOLA and TOLB are the thresholds to determine the effective 00093 * numerical rank of matrix B and a subblock of A. Generally, 00094 * they are set to 00095 * TOLA = MAX(M,N)*norm(A)*MACHEPS, 00096 * TOLB = MAX(P,N)*norm(B)*MACHEPS. 00097 * The size of TOLA and TOLB may affect the size of backward 00098 * errors of the decomposition. 00099 * 00100 * K (output) INTEGER 00101 * L (output) INTEGER 00102 * On exit, K and L specify the dimension of the subblocks 00103 * described in Purpose section. 00104 * K + L = effective numerical rank of (A',B')'. 00105 * 00106 * U (output) COMPLEX array, dimension (LDU,M) 00107 * If JOBU = 'U', U contains the unitary matrix U. 00108 * If JOBU = 'N', U is not referenced. 00109 * 00110 * LDU (input) INTEGER 00111 * The leading dimension of the array U. LDU >= max(1,M) if 00112 * JOBU = 'U'; LDU >= 1 otherwise. 00113 * 00114 * V (output) COMPLEX array, dimension (LDV,P) 00115 * If JOBV = 'V', V contains the unitary matrix V. 00116 * If JOBV = 'N', V is not referenced. 00117 * 00118 * LDV (input) INTEGER 00119 * The leading dimension of the array V. LDV >= max(1,P) if 00120 * JOBV = 'V'; LDV >= 1 otherwise. 00121 * 00122 * Q (output) COMPLEX array, dimension (LDQ,N) 00123 * If JOBQ = 'Q', Q contains the unitary matrix Q. 00124 * If JOBQ = 'N', Q is not referenced. 00125 * 00126 * LDQ (input) INTEGER 00127 * The leading dimension of the array Q. LDQ >= max(1,N) if 00128 * JOBQ = 'Q'; LDQ >= 1 otherwise. 00129 * 00130 * IWORK (workspace) INTEGER array, dimension (N) 00131 * 00132 * RWORK (workspace) REAL array, dimension (2*N) 00133 * 00134 * TAU (workspace) COMPLEX array, dimension (N) 00135 * 00136 * WORK (workspace) COMPLEX array, dimension (max(3*N,M,P)) 00137 * 00138 * INFO (output) INTEGER 00139 * = 0: successful exit 00140 * < 0: if INFO = -i, the i-th argument had an illegal value. 00141 * 00142 * Further Details 00143 * =============== 00144 * 00145 * The subroutine uses LAPACK subroutine CGEQPF for the QR factorization 00146 * with column pivoting to detect the effective numerical rank of the 00147 * a matrix. It may be replaced by a better rank determination strategy. 00148 * 00149 * ===================================================================== 00150 * 00151 * .. Parameters .. 00152 COMPLEX CZERO, CONE 00153 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 00154 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 00155 * .. 00156 * .. Local Scalars .. 00157 LOGICAL FORWRD, WANTQ, WANTU, WANTV 00158 INTEGER I, J 00159 COMPLEX T 00160 * .. 00161 * .. External Functions .. 00162 LOGICAL LSAME 00163 EXTERNAL LSAME 00164 * .. 00165 * .. External Subroutines .. 00166 EXTERNAL CGEQPF, CGEQR2, CGERQ2, CLACPY, CLAPMT, CLASET, 00167 $ CUNG2R, CUNM2R, CUNMR2, XERBLA 00168 * .. 00169 * .. Intrinsic Functions .. 00170 INTRINSIC ABS, AIMAG, MAX, MIN, REAL 00171 * .. 00172 * .. Statement Functions .. 00173 REAL CABS1 00174 * .. 00175 * .. Statement Function definitions .. 00176 CABS1( T ) = ABS( REAL( T ) ) + ABS( AIMAG( T ) ) 00177 * .. 00178 * .. Executable Statements .. 00179 * 00180 * Test the input parameters 00181 * 00182 WANTU = LSAME( JOBU, 'U' ) 00183 WANTV = LSAME( JOBV, 'V' ) 00184 WANTQ = LSAME( JOBQ, 'Q' ) 00185 FORWRD = .TRUE. 00186 * 00187 INFO = 0 00188 IF( .NOT.( WANTU .OR. LSAME( JOBU, 'N' ) ) ) THEN 00189 INFO = -1 00190 ELSE IF( .NOT.( WANTV .OR. LSAME( JOBV, 'N' ) ) ) THEN 00191 INFO = -2 00192 ELSE IF( .NOT.( WANTQ .OR. LSAME( JOBQ, 'N' ) ) ) THEN 00193 INFO = -3 00194 ELSE IF( M.LT.0 ) THEN 00195 INFO = -4 00196 ELSE IF( P.LT.0 ) THEN 00197 INFO = -5 00198 ELSE IF( N.LT.0 ) THEN 00199 INFO = -6 00200 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00201 INFO = -8 00202 ELSE IF( LDB.LT.MAX( 1, P ) ) THEN 00203 INFO = -10 00204 ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN 00205 INFO = -16 00206 ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN 00207 INFO = -18 00208 ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN 00209 INFO = -20 00210 END IF 00211 IF( INFO.NE.0 ) THEN 00212 CALL XERBLA( 'CGGSVP', -INFO ) 00213 RETURN 00214 END IF 00215 * 00216 * QR with column pivoting of B: B*P = V*( S11 S12 ) 00217 * ( 0 0 ) 00218 * 00219 DO 10 I = 1, N 00220 IWORK( I ) = 0 00221 10 CONTINUE 00222 CALL CGEQPF( P, N, B, LDB, IWORK, TAU, WORK, RWORK, INFO ) 00223 * 00224 * Update A := A*P 00225 * 00226 CALL CLAPMT( FORWRD, M, N, A, LDA, IWORK ) 00227 * 00228 * Determine the effective rank of matrix B. 00229 * 00230 L = 0 00231 DO 20 I = 1, MIN( P, N ) 00232 IF( CABS1( B( I, I ) ).GT.TOLB ) 00233 $ L = L + 1 00234 20 CONTINUE 00235 * 00236 IF( WANTV ) THEN 00237 * 00238 * Copy the details of V, and form V. 00239 * 00240 CALL CLASET( 'Full', P, P, CZERO, CZERO, V, LDV ) 00241 IF( P.GT.1 ) 00242 $ CALL CLACPY( 'Lower', P-1, N, B( 2, 1 ), LDB, V( 2, 1 ), 00243 $ LDV ) 00244 CALL CUNG2R( P, P, MIN( P, N ), V, LDV, TAU, WORK, INFO ) 00245 END IF 00246 * 00247 * Clean up B 00248 * 00249 DO 40 J = 1, L - 1 00250 DO 30 I = J + 1, L 00251 B( I, J ) = CZERO 00252 30 CONTINUE 00253 40 CONTINUE 00254 IF( P.GT.L ) 00255 $ CALL CLASET( 'Full', P-L, N, CZERO, CZERO, B( L+1, 1 ), LDB ) 00256 * 00257 IF( WANTQ ) THEN 00258 * 00259 * Set Q = I and Update Q := Q*P 00260 * 00261 CALL CLASET( 'Full', N, N, CZERO, CONE, Q, LDQ ) 00262 CALL CLAPMT( FORWRD, N, N, Q, LDQ, IWORK ) 00263 END IF 00264 * 00265 IF( P.GE.L .AND. N.NE.L ) THEN 00266 * 00267 * RQ factorization of ( S11 S12 ) = ( 0 S12 )*Z 00268 * 00269 CALL CGERQ2( L, N, B, LDB, TAU, WORK, INFO ) 00270 * 00271 * Update A := A*Z' 00272 * 00273 CALL CUNMR2( 'Right', 'Conjugate transpose', M, N, L, B, LDB, 00274 $ TAU, A, LDA, WORK, INFO ) 00275 IF( WANTQ ) THEN 00276 * 00277 * Update Q := Q*Z' 00278 * 00279 CALL CUNMR2( 'Right', 'Conjugate transpose', N, N, L, B, 00280 $ LDB, TAU, Q, LDQ, WORK, INFO ) 00281 END IF 00282 * 00283 * Clean up B 00284 * 00285 CALL CLASET( 'Full', L, N-L, CZERO, CZERO, B, LDB ) 00286 DO 60 J = N - L + 1, N 00287 DO 50 I = J - N + L + 1, L 00288 B( I, J ) = CZERO 00289 50 CONTINUE 00290 60 CONTINUE 00291 * 00292 END IF 00293 * 00294 * Let N-L L 00295 * A = ( A11 A12 ) M, 00296 * 00297 * then the following does the complete QR decomposition of A11: 00298 * 00299 * A11 = U*( 0 T12 )*P1' 00300 * ( 0 0 ) 00301 * 00302 DO 70 I = 1, N - L 00303 IWORK( I ) = 0 00304 70 CONTINUE 00305 CALL CGEQPF( M, N-L, A, LDA, IWORK, TAU, WORK, RWORK, INFO ) 00306 * 00307 * Determine the effective rank of A11 00308 * 00309 K = 0 00310 DO 80 I = 1, MIN( M, N-L ) 00311 IF( CABS1( A( I, I ) ).GT.TOLA ) 00312 $ K = K + 1 00313 80 CONTINUE 00314 * 00315 * Update A12 := U'*A12, where A12 = A( 1:M, N-L+1:N ) 00316 * 00317 CALL CUNM2R( 'Left', 'Conjugate transpose', M, L, MIN( M, N-L ), 00318 $ A, LDA, TAU, A( 1, N-L+1 ), LDA, WORK, INFO ) 00319 * 00320 IF( WANTU ) THEN 00321 * 00322 * Copy the details of U, and form U 00323 * 00324 CALL CLASET( 'Full', M, M, CZERO, CZERO, U, LDU ) 00325 IF( M.GT.1 ) 00326 $ CALL CLACPY( 'Lower', M-1, N-L, A( 2, 1 ), LDA, U( 2, 1 ), 00327 $ LDU ) 00328 CALL CUNG2R( M, M, MIN( M, N-L ), U, LDU, TAU, WORK, INFO ) 00329 END IF 00330 * 00331 IF( WANTQ ) THEN 00332 * 00333 * Update Q( 1:N, 1:N-L ) = Q( 1:N, 1:N-L )*P1 00334 * 00335 CALL CLAPMT( FORWRD, N, N-L, Q, LDQ, IWORK ) 00336 END IF 00337 * 00338 * Clean up A: set the strictly lower triangular part of 00339 * A(1:K, 1:K) = 0, and A( K+1:M, 1:N-L ) = 0. 00340 * 00341 DO 100 J = 1, K - 1 00342 DO 90 I = J + 1, K 00343 A( I, J ) = CZERO 00344 90 CONTINUE 00345 100 CONTINUE 00346 IF( M.GT.K ) 00347 $ CALL CLASET( 'Full', M-K, N-L, CZERO, CZERO, A( K+1, 1 ), LDA ) 00348 * 00349 IF( N-L.GT.K ) THEN 00350 * 00351 * RQ factorization of ( T11 T12 ) = ( 0 T12 )*Z1 00352 * 00353 CALL CGERQ2( K, N-L, A, LDA, TAU, WORK, INFO ) 00354 * 00355 IF( WANTQ ) THEN 00356 * 00357 * Update Q( 1:N,1:N-L ) = Q( 1:N,1:N-L )*Z1' 00358 * 00359 CALL CUNMR2( 'Right', 'Conjugate transpose', N, N-L, K, A, 00360 $ LDA, TAU, Q, LDQ, WORK, INFO ) 00361 END IF 00362 * 00363 * Clean up A 00364 * 00365 CALL CLASET( 'Full', K, N-L-K, CZERO, CZERO, A, LDA ) 00366 DO 120 J = N - L - K + 1, N - L 00367 DO 110 I = J - N + L + K + 1, K 00368 A( I, J ) = CZERO 00369 110 CONTINUE 00370 120 CONTINUE 00371 * 00372 END IF 00373 * 00374 IF( M.GT.K ) THEN 00375 * 00376 * QR factorization of A( K+1:M,N-L+1:N ) 00377 * 00378 CALL CGEQR2( M-K, L, A( K+1, N-L+1 ), LDA, TAU, WORK, INFO ) 00379 * 00380 IF( WANTU ) THEN 00381 * 00382 * Update U(:,K+1:M) := U(:,K+1:M)*U1 00383 * 00384 CALL CUNM2R( 'Right', 'No transpose', M, M-K, MIN( M-K, L ), 00385 $ A( K+1, N-L+1 ), LDA, TAU, U( 1, K+1 ), LDU, 00386 $ WORK, INFO ) 00387 END IF 00388 * 00389 * Clean up 00390 * 00391 DO 140 J = N - L + 1, N 00392 DO 130 I = J - N + K + L + 1, M 00393 A( I, J ) = CZERO 00394 130 CONTINUE 00395 140 CONTINUE 00396 * 00397 END IF 00398 * 00399 RETURN 00400 * 00401 * End of CGGSVP 00402 * 00403 END
1.7.3 
