1*> \brief \b ZGGSVP 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download ZGGSVP + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zggsvp.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zggsvp.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zggsvp.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE ZGGSVP( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, 22* TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ, 23* IWORK, RWORK, TAU, WORK, INFO ) 24* 25* .. Scalar Arguments .. 26* CHARACTER JOBQ, JOBU, JOBV 27* INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P 28* DOUBLE PRECISION TOLA, TOLB 29* .. 30* .. Array Arguments .. 31* INTEGER IWORK( * ) 32* DOUBLE PRECISION RWORK( * ) 33* COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ), 34* $ TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * ) 35* .. 36* 37* 38*> \par Purpose: 39* ============= 40*> 41*> \verbatim 42*> 43*> This routine is deprecated and has been replaced by routine ZGGSVP3. 44*> 45*> ZGGSVP computes unitary matrices U, V and Q such that 46*> 47*> N-K-L K L 48*> U**H*A*Q = K ( 0 A12 A13 ) if M-K-L >= 0; 49*> L ( 0 0 A23 ) 50*> M-K-L ( 0 0 0 ) 51*> 52*> N-K-L K L 53*> = K ( 0 A12 A13 ) if M-K-L < 0; 54*> M-K ( 0 0 A23 ) 55*> 56*> N-K-L K L 57*> V**H*B*Q = L ( 0 0 B13 ) 58*> P-L ( 0 0 0 ) 59*> 60*> where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular 61*> upper triangular; A23 is L-by-L upper triangular if M-K-L >= 0, 62*> otherwise A23 is (M-K)-by-L upper trapezoidal. K+L = the effective 63*> numerical rank of the (M+P)-by-N matrix (A**H,B**H)**H. 64*> 65*> This decomposition is the preprocessing step for computing the 66*> Generalized Singular Value Decomposition (GSVD), see subroutine 67*> ZGGSVD. 68*> \endverbatim 69* 70* Arguments: 71* ========== 72* 73*> \param[in] JOBU 74*> \verbatim 75*> JOBU is CHARACTER*1 76*> = 'U': Unitary matrix U is computed; 77*> = 'N': U is not computed. 78*> \endverbatim 79*> 80*> \param[in] JOBV 81*> \verbatim 82*> JOBV is CHARACTER*1 83*> = 'V': Unitary matrix V is computed; 84*> = 'N': V is not computed. 85*> \endverbatim 86*> 87*> \param[in] JOBQ 88*> \verbatim 89*> JOBQ is CHARACTER*1 90*> = 'Q': Unitary matrix Q is computed; 91*> = 'N': Q is not computed. 92*> \endverbatim 93*> 94*> \param[in] M 95*> \verbatim 96*> M is INTEGER 97*> The number of rows of the matrix A. M >= 0. 98*> \endverbatim 99*> 100*> \param[in] P 101*> \verbatim 102*> P is INTEGER 103*> The number of rows of the matrix B. P >= 0. 104*> \endverbatim 105*> 106*> \param[in] N 107*> \verbatim 108*> N is INTEGER 109*> The number of columns of the matrices A and B. N >= 0. 110*> \endverbatim 111*> 112*> \param[in,out] A 113*> \verbatim 114*> A is COMPLEX*16 array, dimension (LDA,N) 115*> On entry, the M-by-N matrix A. 116*> On exit, A contains the triangular (or trapezoidal) matrix 117*> described in the Purpose section. 118*> \endverbatim 119*> 120*> \param[in] LDA 121*> \verbatim 122*> LDA is INTEGER 123*> The leading dimension of the array A. LDA >= max(1,M). 124*> \endverbatim 125*> 126*> \param[in,out] B 127*> \verbatim 128*> B is COMPLEX*16 array, dimension (LDB,N) 129*> On entry, the P-by-N matrix B. 130*> On exit, B contains the triangular matrix described in 131*> the Purpose section. 132*> \endverbatim 133*> 134*> \param[in] LDB 135*> \verbatim 136*> LDB is INTEGER 137*> The leading dimension of the array B. LDB >= max(1,P). 138*> \endverbatim 139*> 140*> \param[in] TOLA 141*> \verbatim 142*> TOLA is DOUBLE PRECISION 143*> \endverbatim 144*> 145*> \param[in] TOLB 146*> \verbatim 147*> TOLB is DOUBLE PRECISION 148*> 149*> TOLA and TOLB are the thresholds to determine the effective 150*> numerical rank of matrix B and a subblock of A. Generally, 151*> they are set to 152*> TOLA = MAX(M,N)*norm(A)*MAZHEPS, 153*> TOLB = MAX(P,N)*norm(B)*MAZHEPS. 154*> The size of TOLA and TOLB may affect the size of backward 155*> errors of the decomposition. 156*> \endverbatim 157*> 158*> \param[out] K 159*> \verbatim 160*> K is INTEGER 161*> \endverbatim 162*> 163*> \param[out] L 164*> \verbatim 165*> L is INTEGER 166*> 167*> On exit, K and L specify the dimension of the subblocks 168*> described in Purpose section. 169*> K + L = effective numerical rank of (A**H,B**H)**H. 170*> \endverbatim 171*> 172*> \param[out] U 173*> \verbatim 174*> U is COMPLEX*16 array, dimension (LDU,M) 175*> If JOBU = 'U', U contains the unitary matrix U. 176*> If JOBU = 'N', U is not referenced. 177*> \endverbatim 178*> 179*> \param[in] LDU 180*> \verbatim 181*> LDU is INTEGER 182*> The leading dimension of the array U. LDU >= max(1,M) if 183*> JOBU = 'U'; LDU >= 1 otherwise. 184*> \endverbatim 185*> 186*> \param[out] V 187*> \verbatim 188*> V is COMPLEX*16 array, dimension (LDV,P) 189*> If JOBV = 'V', V contains the unitary matrix V. 190*> If JOBV = 'N', V is not referenced. 191*> \endverbatim 192*> 193*> \param[in] LDV 194*> \verbatim 195*> LDV is INTEGER 196*> The leading dimension of the array V. LDV >= max(1,P) if 197*> JOBV = 'V'; LDV >= 1 otherwise. 198*> \endverbatim 199*> 200*> \param[out] Q 201*> \verbatim 202*> Q is COMPLEX*16 array, dimension (LDQ,N) 203*> If JOBQ = 'Q', Q contains the unitary matrix Q. 204*> If JOBQ = 'N', Q is not referenced. 205*> \endverbatim 206*> 207*> \param[in] LDQ 208*> \verbatim 209*> LDQ is INTEGER 210*> The leading dimension of the array Q. LDQ >= max(1,N) if 211*> JOBQ = 'Q'; LDQ >= 1 otherwise. 212*> \endverbatim 213*> 214*> \param[out] IWORK 215*> \verbatim 216*> IWORK is INTEGER array, dimension (N) 217*> \endverbatim 218*> 219*> \param[out] RWORK 220*> \verbatim 221*> RWORK is DOUBLE PRECISION array, dimension (2*N) 222*> \endverbatim 223*> 224*> \param[out] TAU 225*> \verbatim 226*> TAU is COMPLEX*16 array, dimension (N) 227*> \endverbatim 228*> 229*> \param[out] WORK 230*> \verbatim 231*> WORK is COMPLEX*16 array, dimension (max(3*N,M,P)) 232*> \endverbatim 233*> 234*> \param[out] INFO 235*> \verbatim 236*> INFO is INTEGER 237*> = 0: successful exit 238*> < 0: if INFO = -i, the i-th argument had an illegal value. 239*> \endverbatim 240* 241* Authors: 242* ======== 243* 244*> \author Univ. of Tennessee 245*> \author Univ. of California Berkeley 246*> \author Univ. of Colorado Denver 247*> \author NAG Ltd. 248* 249*> \date November 2011 250* 251*> \ingroup complex16OTHERcomputational 252* 253*> \par Further Details: 254* ===================== 255*> 256*> \verbatim 257*> 258*> The subroutine uses LAPACK subroutine ZGEQPF for the QR factorization 259*> with column pivoting to detect the effective numerical rank of the 260*> a matrix. It may be replaced by a better rank determination strategy. 261*> \endverbatim 262*> 263* ===================================================================== 264 SUBROUTINE ZGGSVP( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, 265 $ TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ, 266 $ IWORK, RWORK, TAU, WORK, INFO ) 267* 268* -- LAPACK computational routine (version 3.4.0) -- 269* -- LAPACK is a software package provided by Univ. of Tennessee, -- 270* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 271* November 2011 272* 273* .. Scalar Arguments .. 274 CHARACTER JOBQ, JOBU, JOBV 275 INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P 276 DOUBLE PRECISION TOLA, TOLB 277* .. 278* .. Array Arguments .. 279 INTEGER IWORK( * ) 280 DOUBLE PRECISION RWORK( * ) 281 COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ), 282 $ TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * ) 283* .. 284* 285* ===================================================================== 286* 287* .. Parameters .. 288 COMPLEX*16 CZERO, CONE 289 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 290 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 291* .. 292* .. Local Scalars .. 293 LOGICAL FORWRD, WANTQ, WANTU, WANTV 294 INTEGER I, J 295 COMPLEX*16 T 296* .. 297* .. External Functions .. 298 LOGICAL LSAME 299 EXTERNAL LSAME 300* .. 301* .. External Subroutines .. 302 EXTERNAL XERBLA, ZGEQPF, ZGEQR2, ZGERQ2, ZLACPY, ZLAPMT, 303 $ ZLASET, ZUNG2R, ZUNM2R, ZUNMR2 304* .. 305* .. Intrinsic Functions .. 306 INTRINSIC ABS, DBLE, DIMAG, MAX, MIN 307* .. 308* .. Statement Functions .. 309 DOUBLE PRECISION CABS1 310* .. 311* .. Statement Function definitions .. 312 CABS1( T ) = ABS( DBLE( T ) ) + ABS( DIMAG( T ) ) 313* .. 314* .. Executable Statements .. 315* 316* Test the input parameters 317* 318 WANTU = LSAME( JOBU, 'U' ) 319 WANTV = LSAME( JOBV, 'V' ) 320 WANTQ = LSAME( JOBQ, 'Q' ) 321 FORWRD = .TRUE. 322* 323 INFO = 0 324 IF( .NOT.( WANTU .OR. LSAME( JOBU, 'N' ) ) ) THEN 325 INFO = -1 326 ELSE IF( .NOT.( WANTV .OR. LSAME( JOBV, 'N' ) ) ) THEN 327 INFO = -2 328 ELSE IF( .NOT.( WANTQ .OR. LSAME( JOBQ, 'N' ) ) ) THEN 329 INFO = -3 330 ELSE IF( M.LT.0 ) THEN 331 INFO = -4 332 ELSE IF( P.LT.0 ) THEN 333 INFO = -5 334 ELSE IF( N.LT.0 ) THEN 335 INFO = -6 336 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 337 INFO = -8 338 ELSE IF( LDB.LT.MAX( 1, P ) ) THEN 339 INFO = -10 340 ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN 341 INFO = -16 342 ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN 343 INFO = -18 344 ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN 345 INFO = -20 346 END IF 347 IF( INFO.NE.0 ) THEN 348 CALL XERBLA( 'ZGGSVP', -INFO ) 349 RETURN 350 END IF 351* 352* QR with column pivoting of B: B*P = V*( S11 S12 ) 353* ( 0 0 ) 354* 355 DO 10 I = 1, N 356 IWORK( I ) = 0 357 10 CONTINUE 358 CALL ZGEQPF( P, N, B, LDB, IWORK, TAU, WORK, RWORK, INFO ) 359* 360* Update A := A*P 361* 362 CALL ZLAPMT( FORWRD, M, N, A, LDA, IWORK ) 363* 364* Determine the effective rank of matrix B. 365* 366 L = 0 367 DO 20 I = 1, MIN( P, N ) 368 IF( CABS1( B( I, I ) ).GT.TOLB ) 369 $ L = L + 1 370 20 CONTINUE 371* 372 IF( WANTV ) THEN 373* 374* Copy the details of V, and form V. 375* 376 CALL ZLASET( 'Full', P, P, CZERO, CZERO, V, LDV ) 377 IF( P.GT.1 ) 378 $ CALL ZLACPY( 'Lower', P-1, N, B( 2, 1 ), LDB, V( 2, 1 ), 379 $ LDV ) 380 CALL ZUNG2R( P, P, MIN( P, N ), V, LDV, TAU, WORK, INFO ) 381 END IF 382* 383* Clean up B 384* 385 DO 40 J = 1, L - 1 386 DO 30 I = J + 1, L 387 B( I, J ) = CZERO 388 30 CONTINUE 389 40 CONTINUE 390 IF( P.GT.L ) 391 $ CALL ZLASET( 'Full', P-L, N, CZERO, CZERO, B( L+1, 1 ), LDB ) 392* 393 IF( WANTQ ) THEN 394* 395* Set Q = I and Update Q := Q*P 396* 397 CALL ZLASET( 'Full', N, N, CZERO, CONE, Q, LDQ ) 398 CALL ZLAPMT( FORWRD, N, N, Q, LDQ, IWORK ) 399 END IF 400* 401 IF( P.GE.L .AND. N.NE.L ) THEN 402* 403* RQ factorization of ( S11 S12 ) = ( 0 S12 )*Z 404* 405 CALL ZGERQ2( L, N, B, LDB, TAU, WORK, INFO ) 406* 407* Update A := A*Z**H 408* 409 CALL ZUNMR2( 'Right', 'Conjugate transpose', M, N, L, B, LDB, 410 $ TAU, A, LDA, WORK, INFO ) 411 IF( WANTQ ) THEN 412* 413* Update Q := Q*Z**H 414* 415 CALL ZUNMR2( 'Right', 'Conjugate transpose', N, N, L, B, 416 $ LDB, TAU, Q, LDQ, WORK, INFO ) 417 END IF 418* 419* Clean up B 420* 421 CALL ZLASET( 'Full', L, N-L, CZERO, CZERO, B, LDB ) 422 DO 60 J = N - L + 1, N 423 DO 50 I = J - N + L + 1, L 424 B( I, J ) = CZERO 425 50 CONTINUE 426 60 CONTINUE 427* 428 END IF 429* 430* Let N-L L 431* A = ( A11 A12 ) M, 432* 433* then the following does the complete QR decomposition of A11: 434* 435* A11 = U*( 0 T12 )*P1**H 436* ( 0 0 ) 437* 438 DO 70 I = 1, N - L 439 IWORK( I ) = 0 440 70 CONTINUE 441 CALL ZGEQPF( M, N-L, A, LDA, IWORK, TAU, WORK, RWORK, INFO ) 442* 443* Determine the effective rank of A11 444* 445 K = 0 446 DO 80 I = 1, MIN( M, N-L ) 447 IF( CABS1( A( I, I ) ).GT.TOLA ) 448 $ K = K + 1 449 80 CONTINUE 450* 451* Update A12 := U**H*A12, where A12 = A( 1:M, N-L+1:N ) 452* 453 CALL ZUNM2R( 'Left', 'Conjugate transpose', M, L, MIN( M, N-L ), 454 $ A, LDA, TAU, A( 1, N-L+1 ), LDA, WORK, INFO ) 455* 456 IF( WANTU ) THEN 457* 458* Copy the details of U, and form U 459* 460 CALL ZLASET( 'Full', M, M, CZERO, CZERO, U, LDU ) 461 IF( M.GT.1 ) 462 $ CALL ZLACPY( 'Lower', M-1, N-L, A( 2, 1 ), LDA, U( 2, 1 ), 463 $ LDU ) 464 CALL ZUNG2R( M, M, MIN( M, N-L ), U, LDU, TAU, WORK, INFO ) 465 END IF 466* 467 IF( WANTQ ) THEN 468* 469* Update Q( 1:N, 1:N-L ) = Q( 1:N, 1:N-L )*P1 470* 471 CALL ZLAPMT( FORWRD, N, N-L, Q, LDQ, IWORK ) 472 END IF 473* 474* Clean up A: set the strictly lower triangular part of 475* A(1:K, 1:K) = 0, and A( K+1:M, 1:N-L ) = 0. 476* 477 DO 100 J = 1, K - 1 478 DO 90 I = J + 1, K 479 A( I, J ) = CZERO 480 90 CONTINUE 481 100 CONTINUE 482 IF( M.GT.K ) 483 $ CALL ZLASET( 'Full', M-K, N-L, CZERO, CZERO, A( K+1, 1 ), LDA ) 484* 485 IF( N-L.GT.K ) THEN 486* 487* RQ factorization of ( T11 T12 ) = ( 0 T12 )*Z1 488* 489 CALL ZGERQ2( K, N-L, A, LDA, TAU, WORK, INFO ) 490* 491 IF( WANTQ ) THEN 492* 493* Update Q( 1:N,1:N-L ) = Q( 1:N,1:N-L )*Z1**H 494* 495 CALL ZUNMR2( 'Right', 'Conjugate transpose', N, N-L, K, A, 496 $ LDA, TAU, Q, LDQ, WORK, INFO ) 497 END IF 498* 499* Clean up A 500* 501 CALL ZLASET( 'Full', K, N-L-K, CZERO, CZERO, A, LDA ) 502 DO 120 J = N - L - K + 1, N - L 503 DO 110 I = J - N + L + K + 1, K 504 A( I, J ) = CZERO 505 110 CONTINUE 506 120 CONTINUE 507* 508 END IF 509* 510 IF( M.GT.K ) THEN 511* 512* QR factorization of A( K+1:M,N-L+1:N ) 513* 514 CALL ZGEQR2( M-K, L, A( K+1, N-L+1 ), LDA, TAU, WORK, INFO ) 515* 516 IF( WANTU ) THEN 517* 518* Update U(:,K+1:M) := U(:,K+1:M)*U1 519* 520 CALL ZUNM2R( 'Right', 'No transpose', M, M-K, MIN( M-K, L ), 521 $ A( K+1, N-L+1 ), LDA, TAU, U( 1, K+1 ), LDU, 522 $ WORK, INFO ) 523 END IF 524* 525* Clean up 526* 527 DO 140 J = N - L + 1, N 528 DO 130 I = J - N + K + L + 1, M 529 A( I, J ) = CZERO 530 130 CONTINUE 531 140 CONTINUE 532* 533 END IF 534* 535 RETURN 536* 537* End of ZGGSVP 538* 539 END 540