![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZGGSVP 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZGGSVP + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zggsvp.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zggsvp.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zggsvp.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZGGSVP( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, 00022 * TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ, 00023 * IWORK, RWORK, TAU, WORK, INFO ) 00024 * 00025 * .. Scalar Arguments .. 00026 * CHARACTER JOBQ, JOBU, JOBV 00027 * INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P 00028 * DOUBLE PRECISION TOLA, TOLB 00029 * .. 00030 * .. Array Arguments .. 00031 * INTEGER IWORK( * ) 00032 * DOUBLE PRECISION RWORK( * ) 00033 * COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ), 00034 * $ TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * ) 00035 * .. 00036 * 00037 * 00038 *> \par Purpose: 00039 * ============= 00040 *> 00041 *> \verbatim 00042 *> 00043 *> ZGGSVP computes unitary matrices U, V and Q such that 00044 *> 00045 *> N-K-L K L 00046 *> U**H*A*Q = K ( 0 A12 A13 ) if M-K-L >= 0; 00047 *> L ( 0 0 A23 ) 00048 *> M-K-L ( 0 0 0 ) 00049 *> 00050 *> N-K-L K L 00051 *> = K ( 0 A12 A13 ) if M-K-L < 0; 00052 *> M-K ( 0 0 A23 ) 00053 *> 00054 *> N-K-L K L 00055 *> V**H*B*Q = L ( 0 0 B13 ) 00056 *> P-L ( 0 0 0 ) 00057 *> 00058 *> where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular 00059 *> upper triangular; A23 is L-by-L upper triangular if M-K-L >= 0, 00060 *> otherwise A23 is (M-K)-by-L upper trapezoidal. K+L = the effective 00061 *> numerical rank of the (M+P)-by-N matrix (A**H,B**H)**H. 00062 *> 00063 *> This decomposition is the preprocessing step for computing the 00064 *> Generalized Singular Value Decomposition (GSVD), see subroutine 00065 *> ZGGSVD. 00066 *> \endverbatim 00067 * 00068 * Arguments: 00069 * ========== 00070 * 00071 *> \param[in] JOBU 00072 *> \verbatim 00073 *> JOBU is CHARACTER*1 00074 *> = 'U': Unitary matrix U is computed; 00075 *> = 'N': U is not computed. 00076 *> \endverbatim 00077 *> 00078 *> \param[in] JOBV 00079 *> \verbatim 00080 *> JOBV is CHARACTER*1 00081 *> = 'V': Unitary matrix V is computed; 00082 *> = 'N': V is not computed. 00083 *> \endverbatim 00084 *> 00085 *> \param[in] JOBQ 00086 *> \verbatim 00087 *> JOBQ is CHARACTER*1 00088 *> = 'Q': Unitary matrix Q is computed; 00089 *> = 'N': Q is not computed. 00090 *> \endverbatim 00091 *> 00092 *> \param[in] M 00093 *> \verbatim 00094 *> M is INTEGER 00095 *> The number of rows of the matrix A. M >= 0. 00096 *> \endverbatim 00097 *> 00098 *> \param[in] P 00099 *> \verbatim 00100 *> P is INTEGER 00101 *> The number of rows of the matrix B. P >= 0. 00102 *> \endverbatim 00103 *> 00104 *> \param[in] N 00105 *> \verbatim 00106 *> N is INTEGER 00107 *> The number of columns of the matrices A and B. N >= 0. 00108 *> \endverbatim 00109 *> 00110 *> \param[in,out] A 00111 *> \verbatim 00112 *> A is COMPLEX*16 array, dimension (LDA,N) 00113 *> On entry, the M-by-N matrix A. 00114 *> On exit, A contains the triangular (or trapezoidal) matrix 00115 *> described in the Purpose section. 00116 *> \endverbatim 00117 *> 00118 *> \param[in] LDA 00119 *> \verbatim 00120 *> LDA is INTEGER 00121 *> The leading dimension of the array A. LDA >= max(1,M). 00122 *> \endverbatim 00123 *> 00124 *> \param[in,out] B 00125 *> \verbatim 00126 *> B is COMPLEX*16 array, dimension (LDB,N) 00127 *> On entry, the P-by-N matrix B. 00128 *> On exit, B contains the triangular matrix described in 00129 *> the Purpose section. 00130 *> \endverbatim 00131 *> 00132 *> \param[in] LDB 00133 *> \verbatim 00134 *> LDB is INTEGER 00135 *> The leading dimension of the array B. LDB >= max(1,P). 00136 *> \endverbatim 00137 *> 00138 *> \param[in] TOLA 00139 *> \verbatim 00140 *> TOLA is DOUBLE PRECISION 00141 *> \endverbatim 00142 *> 00143 *> \param[in] TOLB 00144 *> \verbatim 00145 *> TOLB is DOUBLE PRECISION 00146 *> 00147 *> TOLA and TOLB are the thresholds to determine the effective 00148 *> numerical rank of matrix B and a subblock of A. Generally, 00149 *> they are set to 00150 *> TOLA = MAX(M,N)*norm(A)*MAZHEPS, 00151 *> TOLB = MAX(P,N)*norm(B)*MAZHEPS. 00152 *> The size of TOLA and TOLB may affect the size of backward 00153 *> errors of the decomposition. 00154 *> \endverbatim 00155 *> 00156 *> \param[out] K 00157 *> \verbatim 00158 *> K is INTEGER 00159 *> \endverbatim 00160 *> 00161 *> \param[out] L 00162 *> \verbatim 00163 *> L is INTEGER 00164 *> 00165 *> On exit, K and L specify the dimension of the subblocks 00166 *> described in Purpose section. 00167 *> K + L = effective numerical rank of (A**H,B**H)**H. 00168 *> \endverbatim 00169 *> 00170 *> \param[out] U 00171 *> \verbatim 00172 *> U is COMPLEX*16 array, dimension (LDU,M) 00173 *> If JOBU = 'U', U contains the unitary matrix U. 00174 *> If JOBU = 'N', U is not referenced. 00175 *> \endverbatim 00176 *> 00177 *> \param[in] LDU 00178 *> \verbatim 00179 *> LDU is INTEGER 00180 *> The leading dimension of the array U. LDU >= max(1,M) if 00181 *> JOBU = 'U'; LDU >= 1 otherwise. 00182 *> \endverbatim 00183 *> 00184 *> \param[out] V 00185 *> \verbatim 00186 *> V is COMPLEX*16 array, dimension (LDV,P) 00187 *> If JOBV = 'V', V contains the unitary matrix V. 00188 *> If JOBV = 'N', V is not referenced. 00189 *> \endverbatim 00190 *> 00191 *> \param[in] LDV 00192 *> \verbatim 00193 *> LDV is INTEGER 00194 *> The leading dimension of the array V. LDV >= max(1,P) if 00195 *> JOBV = 'V'; LDV >= 1 otherwise. 00196 *> \endverbatim 00197 *> 00198 *> \param[out] Q 00199 *> \verbatim 00200 *> Q is COMPLEX*16 array, dimension (LDQ,N) 00201 *> If JOBQ = 'Q', Q contains the unitary matrix Q. 00202 *> If JOBQ = 'N', Q is not referenced. 00203 *> \endverbatim 00204 *> 00205 *> \param[in] LDQ 00206 *> \verbatim 00207 *> LDQ is INTEGER 00208 *> The leading dimension of the array Q. LDQ >= max(1,N) if 00209 *> JOBQ = 'Q'; LDQ >= 1 otherwise. 00210 *> \endverbatim 00211 *> 00212 *> \param[out] IWORK 00213 *> \verbatim 00214 *> IWORK is INTEGER array, dimension (N) 00215 *> \endverbatim 00216 *> 00217 *> \param[out] RWORK 00218 *> \verbatim 00219 *> RWORK is DOUBLE PRECISION array, dimension (2*N) 00220 *> \endverbatim 00221 *> 00222 *> \param[out] TAU 00223 *> \verbatim 00224 *> TAU is COMPLEX*16 array, dimension (N) 00225 *> \endverbatim 00226 *> 00227 *> \param[out] WORK 00228 *> \verbatim 00229 *> WORK is COMPLEX*16 array, dimension (max(3*N,M,P)) 00230 *> \endverbatim 00231 *> 00232 *> \param[out] INFO 00233 *> \verbatim 00234 *> INFO is INTEGER 00235 *> = 0: successful exit 00236 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00237 *> \endverbatim 00238 * 00239 * Authors: 00240 * ======== 00241 * 00242 *> \author Univ. of Tennessee 00243 *> \author Univ. of California Berkeley 00244 *> \author Univ. of Colorado Denver 00245 *> \author NAG Ltd. 00246 * 00247 *> \date November 2011 00248 * 00249 *> \ingroup complex16OTHERcomputational 00250 * 00251 *> \par Further Details: 00252 * ===================== 00253 *> 00254 *> \verbatim 00255 *> 00256 *> The subroutine uses LAPACK subroutine ZGEQPF for the QR factorization 00257 *> with column pivoting to detect the effective numerical rank of the 00258 *> a matrix. It may be replaced by a better rank determination strategy. 00259 *> \endverbatim 00260 *> 00261 * ===================================================================== 00262 SUBROUTINE ZGGSVP( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, 00263 $ TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ, 00264 $ IWORK, RWORK, TAU, WORK, INFO ) 00265 * 00266 * -- LAPACK computational routine (version 3.4.0) -- 00267 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00268 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00269 * November 2011 00270 * 00271 * .. Scalar Arguments .. 00272 CHARACTER JOBQ, JOBU, JOBV 00273 INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P 00274 DOUBLE PRECISION TOLA, TOLB 00275 * .. 00276 * .. Array Arguments .. 00277 INTEGER IWORK( * ) 00278 DOUBLE PRECISION RWORK( * ) 00279 COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ), 00280 $ TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * ) 00281 * .. 00282 * 00283 * ===================================================================== 00284 * 00285 * .. Parameters .. 00286 COMPLEX*16 CZERO, CONE 00287 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 00288 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 00289 * .. 00290 * .. Local Scalars .. 00291 LOGICAL FORWRD, WANTQ, WANTU, WANTV 00292 INTEGER I, J 00293 COMPLEX*16 T 00294 * .. 00295 * .. External Functions .. 00296 LOGICAL LSAME 00297 EXTERNAL LSAME 00298 * .. 00299 * .. External Subroutines .. 00300 EXTERNAL XERBLA, ZGEQPF, ZGEQR2, ZGERQ2, ZLACPY, ZLAPMT, 00301 $ ZLASET, ZUNG2R, ZUNM2R, ZUNMR2 00302 * .. 00303 * .. Intrinsic Functions .. 00304 INTRINSIC ABS, DBLE, DIMAG, MAX, MIN 00305 * .. 00306 * .. Statement Functions .. 00307 DOUBLE PRECISION CABS1 00308 * .. 00309 * .. Statement Function definitions .. 00310 CABS1( T ) = ABS( DBLE( T ) ) + ABS( DIMAG( T ) ) 00311 * .. 00312 * .. Executable Statements .. 00313 * 00314 * Test the input parameters 00315 * 00316 WANTU = LSAME( JOBU, 'U' ) 00317 WANTV = LSAME( JOBV, 'V' ) 00318 WANTQ = LSAME( JOBQ, 'Q' ) 00319 FORWRD = .TRUE. 00320 * 00321 INFO = 0 00322 IF( .NOT.( WANTU .OR. LSAME( JOBU, 'N' ) ) ) THEN 00323 INFO = -1 00324 ELSE IF( .NOT.( WANTV .OR. LSAME( JOBV, 'N' ) ) ) THEN 00325 INFO = -2 00326 ELSE IF( .NOT.( WANTQ .OR. LSAME( JOBQ, 'N' ) ) ) THEN 00327 INFO = -3 00328 ELSE IF( M.LT.0 ) THEN 00329 INFO = -4 00330 ELSE IF( P.LT.0 ) THEN 00331 INFO = -5 00332 ELSE IF( N.LT.0 ) THEN 00333 INFO = -6 00334 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00335 INFO = -8 00336 ELSE IF( LDB.LT.MAX( 1, P ) ) THEN 00337 INFO = -10 00338 ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN 00339 INFO = -16 00340 ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN 00341 INFO = -18 00342 ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN 00343 INFO = -20 00344 END IF 00345 IF( INFO.NE.0 ) THEN 00346 CALL XERBLA( 'ZGGSVP', -INFO ) 00347 RETURN 00348 END IF 00349 * 00350 * QR with column pivoting of B: B*P = V*( S11 S12 ) 00351 * ( 0 0 ) 00352 * 00353 DO 10 I = 1, N 00354 IWORK( I ) = 0 00355 10 CONTINUE 00356 CALL ZGEQPF( P, N, B, LDB, IWORK, TAU, WORK, RWORK, INFO ) 00357 * 00358 * Update A := A*P 00359 * 00360 CALL ZLAPMT( FORWRD, M, N, A, LDA, IWORK ) 00361 * 00362 * Determine the effective rank of matrix B. 00363 * 00364 L = 0 00365 DO 20 I = 1, MIN( P, N ) 00366 IF( CABS1( B( I, I ) ).GT.TOLB ) 00367 $ L = L + 1 00368 20 CONTINUE 00369 * 00370 IF( WANTV ) THEN 00371 * 00372 * Copy the details of V, and form V. 00373 * 00374 CALL ZLASET( 'Full', P, P, CZERO, CZERO, V, LDV ) 00375 IF( P.GT.1 ) 00376 $ CALL ZLACPY( 'Lower', P-1, N, B( 2, 1 ), LDB, V( 2, 1 ), 00377 $ LDV ) 00378 CALL ZUNG2R( P, P, MIN( P, N ), V, LDV, TAU, WORK, INFO ) 00379 END IF 00380 * 00381 * Clean up B 00382 * 00383 DO 40 J = 1, L - 1 00384 DO 30 I = J + 1, L 00385 B( I, J ) = CZERO 00386 30 CONTINUE 00387 40 CONTINUE 00388 IF( P.GT.L ) 00389 $ CALL ZLASET( 'Full', P-L, N, CZERO, CZERO, B( L+1, 1 ), LDB ) 00390 * 00391 IF( WANTQ ) THEN 00392 * 00393 * Set Q = I and Update Q := Q*P 00394 * 00395 CALL ZLASET( 'Full', N, N, CZERO, CONE, Q, LDQ ) 00396 CALL ZLAPMT( FORWRD, N, N, Q, LDQ, IWORK ) 00397 END IF 00398 * 00399 IF( P.GE.L .AND. N.NE.L ) THEN 00400 * 00401 * RQ factorization of ( S11 S12 ) = ( 0 S12 )*Z 00402 * 00403 CALL ZGERQ2( L, N, B, LDB, TAU, WORK, INFO ) 00404 * 00405 * Update A := A*Z**H 00406 * 00407 CALL ZUNMR2( 'Right', 'Conjugate transpose', M, N, L, B, LDB, 00408 $ TAU, A, LDA, WORK, INFO ) 00409 IF( WANTQ ) THEN 00410 * 00411 * Update Q := Q*Z**H 00412 * 00413 CALL ZUNMR2( 'Right', 'Conjugate transpose', N, N, L, B, 00414 $ LDB, TAU, Q, LDQ, WORK, INFO ) 00415 END IF 00416 * 00417 * Clean up B 00418 * 00419 CALL ZLASET( 'Full', L, N-L, CZERO, CZERO, B, LDB ) 00420 DO 60 J = N - L + 1, N 00421 DO 50 I = J - N + L + 1, L 00422 B( I, J ) = CZERO 00423 50 CONTINUE 00424 60 CONTINUE 00425 * 00426 END IF 00427 * 00428 * Let N-L L 00429 * A = ( A11 A12 ) M, 00430 * 00431 * then the following does the complete QR decomposition of A11: 00432 * 00433 * A11 = U*( 0 T12 )*P1**H 00434 * ( 0 0 ) 00435 * 00436 DO 70 I = 1, N - L 00437 IWORK( I ) = 0 00438 70 CONTINUE 00439 CALL ZGEQPF( M, N-L, A, LDA, IWORK, TAU, WORK, RWORK, INFO ) 00440 * 00441 * Determine the effective rank of A11 00442 * 00443 K = 0 00444 DO 80 I = 1, MIN( M, N-L ) 00445 IF( CABS1( A( I, I ) ).GT.TOLA ) 00446 $ K = K + 1 00447 80 CONTINUE 00448 * 00449 * Update A12 := U**H*A12, where A12 = A( 1:M, N-L+1:N ) 00450 * 00451 CALL ZUNM2R( 'Left', 'Conjugate transpose', M, L, MIN( M, N-L ), 00452 $ A, LDA, TAU, A( 1, N-L+1 ), LDA, WORK, INFO ) 00453 * 00454 IF( WANTU ) THEN 00455 * 00456 * Copy the details of U, and form U 00457 * 00458 CALL ZLASET( 'Full', M, M, CZERO, CZERO, U, LDU ) 00459 IF( M.GT.1 ) 00460 $ CALL ZLACPY( 'Lower', M-1, N-L, A( 2, 1 ), LDA, U( 2, 1 ), 00461 $ LDU ) 00462 CALL ZUNG2R( M, M, MIN( M, N-L ), U, LDU, TAU, WORK, INFO ) 00463 END IF 00464 * 00465 IF( WANTQ ) THEN 00466 * 00467 * Update Q( 1:N, 1:N-L ) = Q( 1:N, 1:N-L )*P1 00468 * 00469 CALL ZLAPMT( FORWRD, N, N-L, Q, LDQ, IWORK ) 00470 END IF 00471 * 00472 * Clean up A: set the strictly lower triangular part of 00473 * A(1:K, 1:K) = 0, and A( K+1:M, 1:N-L ) = 0. 00474 * 00475 DO 100 J = 1, K - 1 00476 DO 90 I = J + 1, K 00477 A( I, J ) = CZERO 00478 90 CONTINUE 00479 100 CONTINUE 00480 IF( M.GT.K ) 00481 $ CALL ZLASET( 'Full', M-K, N-L, CZERO, CZERO, A( K+1, 1 ), LDA ) 00482 * 00483 IF( N-L.GT.K ) THEN 00484 * 00485 * RQ factorization of ( T11 T12 ) = ( 0 T12 )*Z1 00486 * 00487 CALL ZGERQ2( K, N-L, A, LDA, TAU, WORK, INFO ) 00488 * 00489 IF( WANTQ ) THEN 00490 * 00491 * Update Q( 1:N,1:N-L ) = Q( 1:N,1:N-L )*Z1**H 00492 * 00493 CALL ZUNMR2( 'Right', 'Conjugate transpose', N, N-L, K, A, 00494 $ LDA, TAU, Q, LDQ, WORK, INFO ) 00495 END IF 00496 * 00497 * Clean up A 00498 * 00499 CALL ZLASET( 'Full', K, N-L-K, CZERO, CZERO, A, LDA ) 00500 DO 120 J = N - L - K + 1, N - L 00501 DO 110 I = J - N + L + K + 1, K 00502 A( I, J ) = CZERO 00503 110 CONTINUE 00504 120 CONTINUE 00505 * 00506 END IF 00507 * 00508 IF( M.GT.K ) THEN 00509 * 00510 * QR factorization of A( K+1:M,N-L+1:N ) 00511 * 00512 CALL ZGEQR2( M-K, L, A( K+1, N-L+1 ), LDA, TAU, WORK, INFO ) 00513 * 00514 IF( WANTU ) THEN 00515 * 00516 * Update U(:,K+1:M) := U(:,K+1:M)*U1 00517 * 00518 CALL ZUNM2R( 'Right', 'No transpose', M, M-K, MIN( M-K, L ), 00519 $ A( K+1, N-L+1 ), LDA, TAU, U( 1, K+1 ), LDU, 00520 $ WORK, INFO ) 00521 END IF 00522 * 00523 * Clean up 00524 * 00525 DO 140 J = N - L + 1, N 00526 DO 130 I = J - N + K + L + 1, M 00527 A( I, J ) = CZERO 00528 130 CONTINUE 00529 140 CONTINUE 00530 * 00531 END IF 00532 * 00533 RETURN 00534 * 00535 * End of ZGGSVP 00536 * 00537 END