1*> \brief \b CUNCSD 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download CUNCSD + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cuncsd.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cuncsd.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cuncsd.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* RECURSIVE SUBROUTINE CUNCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, 22* SIGNS, M, P, Q, X11, LDX11, X12, 23* LDX12, X21, LDX21, X22, LDX22, THETA, 24* U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, 25* LDV2T, WORK, LWORK, RWORK, LRWORK, 26* IWORK, INFO ) 27* 28* .. Scalar Arguments .. 29* CHARACTER JOBU1, JOBU2, JOBV1T, JOBV2T, SIGNS, TRANS 30* INTEGER INFO, LDU1, LDU2, LDV1T, LDV2T, LDX11, LDX12, 31* $ LDX21, LDX22, LRWORK, LWORK, M, P, Q 32* .. 33* .. Array Arguments .. 34* INTEGER IWORK( * ) 35* REAL THETA( * ) 36* REAL RWORK( * ) 37* COMPLEX U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ), 38* $ V2T( LDV2T, * ), WORK( * ), X11( LDX11, * ), 39* $ X12( LDX12, * ), X21( LDX21, * ), X22( LDX22, 40* $ * ) 41* .. 42* 43* 44*> \par Purpose: 45* ============= 46*> 47*> \verbatim 48*> 49*> CUNCSD computes the CS decomposition of an M-by-M partitioned 50*> unitary matrix X: 51*> 52*> [ I 0 0 | 0 0 0 ] 53*> [ 0 C 0 | 0 -S 0 ] 54*> [ X11 | X12 ] [ U1 | ] [ 0 0 0 | 0 0 -I ] [ V1 | ]**H 55*> X = [-----------] = [---------] [---------------------] [---------] . 56*> [ X21 | X22 ] [ | U2 ] [ 0 0 0 | I 0 0 ] [ | V2 ] 57*> [ 0 S 0 | 0 C 0 ] 58*> [ 0 0 I | 0 0 0 ] 59*> 60*> X11 is P-by-Q. The unitary matrices U1, U2, V1, and V2 are P-by-P, 61*> (M-P)-by-(M-P), Q-by-Q, and (M-Q)-by-(M-Q), respectively. C and S are 62*> R-by-R nonnegative diagonal matrices satisfying C^2 + S^2 = I, in 63*> which R = MIN(P,M-P,Q,M-Q). 64*> \endverbatim 65* 66* Arguments: 67* ========== 68* 69*> \param[in] JOBU1 70*> \verbatim 71*> JOBU1 is CHARACTER 72*> = 'Y': U1 is computed; 73*> otherwise: U1 is not computed. 74*> \endverbatim 75*> 76*> \param[in] JOBU2 77*> \verbatim 78*> JOBU2 is CHARACTER 79*> = 'Y': U2 is computed; 80*> otherwise: U2 is not computed. 81*> \endverbatim 82*> 83*> \param[in] JOBV1T 84*> \verbatim 85*> JOBV1T is CHARACTER 86*> = 'Y': V1T is computed; 87*> otherwise: V1T is not computed. 88*> \endverbatim 89*> 90*> \param[in] JOBV2T 91*> \verbatim 92*> JOBV2T is CHARACTER 93*> = 'Y': V2T is computed; 94*> otherwise: V2T is not computed. 95*> \endverbatim 96*> 97*> \param[in] TRANS 98*> \verbatim 99*> TRANS is CHARACTER 100*> = 'T': X, U1, U2, V1T, and V2T are stored in row-major 101*> order; 102*> otherwise: X, U1, U2, V1T, and V2T are stored in column- 103*> major order. 104*> \endverbatim 105*> 106*> \param[in] SIGNS 107*> \verbatim 108*> SIGNS is CHARACTER 109*> = 'O': The lower-left block is made nonpositive (the 110*> "other" convention); 111*> otherwise: The upper-right block is made nonpositive (the 112*> "default" convention). 113*> \endverbatim 114*> 115*> \param[in] M 116*> \verbatim 117*> M is INTEGER 118*> The number of rows and columns in X. 119*> \endverbatim 120*> 121*> \param[in] P 122*> \verbatim 123*> P is INTEGER 124*> The number of rows in X11 and X12. 0 <= P <= M. 125*> \endverbatim 126*> 127*> \param[in] Q 128*> \verbatim 129*> Q is INTEGER 130*> The number of columns in X11 and X21. 0 <= Q <= M. 131*> \endverbatim 132*> 133*> \param[in,out] X11 134*> \verbatim 135*> X11 is COMPLEX array, dimension (LDX11,Q) 136*> On entry, part of the unitary matrix whose CSD is desired. 137*> \endverbatim 138*> 139*> \param[in] LDX11 140*> \verbatim 141*> LDX11 is INTEGER 142*> The leading dimension of X11. LDX11 >= MAX(1,P). 143*> \endverbatim 144*> 145*> \param[in,out] X12 146*> \verbatim 147*> X12 is COMPLEX array, dimension (LDX12,M-Q) 148*> On entry, part of the unitary matrix whose CSD is desired. 149*> \endverbatim 150*> 151*> \param[in] LDX12 152*> \verbatim 153*> LDX12 is INTEGER 154*> The leading dimension of X12. LDX12 >= MAX(1,P). 155*> \endverbatim 156*> 157*> \param[in,out] X21 158*> \verbatim 159*> X21 is COMPLEX array, dimension (LDX21,Q) 160*> On entry, part of the unitary matrix whose CSD is desired. 161*> \endverbatim 162*> 163*> \param[in] LDX21 164*> \verbatim 165*> LDX21 is INTEGER 166*> The leading dimension of X11. LDX21 >= MAX(1,M-P). 167*> \endverbatim 168*> 169*> \param[in,out] X22 170*> \verbatim 171*> X22 is COMPLEX array, dimension (LDX22,M-Q) 172*> On entry, part of the unitary matrix whose CSD is desired. 173*> \endverbatim 174*> 175*> \param[in] LDX22 176*> \verbatim 177*> LDX22 is INTEGER 178*> The leading dimension of X11. LDX22 >= MAX(1,M-P). 179*> \endverbatim 180*> 181*> \param[out] THETA 182*> \verbatim 183*> THETA is REAL array, dimension (R), in which R = 184*> MIN(P,M-P,Q,M-Q). 185*> C = DIAG( COS(THETA(1)), ... , COS(THETA(R)) ) and 186*> S = DIAG( SIN(THETA(1)), ... , SIN(THETA(R)) ). 187*> \endverbatim 188*> 189*> \param[out] U1 190*> \verbatim 191*> U1 is COMPLEX array, dimension (LDU1,P) 192*> If JOBU1 = 'Y', U1 contains the P-by-P unitary matrix U1. 193*> \endverbatim 194*> 195*> \param[in] LDU1 196*> \verbatim 197*> LDU1 is INTEGER 198*> The leading dimension of U1. If JOBU1 = 'Y', LDU1 >= 199*> MAX(1,P). 200*> \endverbatim 201*> 202*> \param[out] U2 203*> \verbatim 204*> U2 is COMPLEX array, dimension (LDU2,M-P) 205*> If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) unitary 206*> matrix U2. 207*> \endverbatim 208*> 209*> \param[in] LDU2 210*> \verbatim 211*> LDU2 is INTEGER 212*> The leading dimension of U2. If JOBU2 = 'Y', LDU2 >= 213*> MAX(1,M-P). 214*> \endverbatim 215*> 216*> \param[out] V1T 217*> \verbatim 218*> V1T is COMPLEX array, dimension (LDV1T,Q) 219*> If JOBV1T = 'Y', V1T contains the Q-by-Q matrix unitary 220*> matrix V1**H. 221*> \endverbatim 222*> 223*> \param[in] LDV1T 224*> \verbatim 225*> LDV1T is INTEGER 226*> The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >= 227*> MAX(1,Q). 228*> \endverbatim 229*> 230*> \param[out] V2T 231*> \verbatim 232*> V2T is COMPLEX array, dimension (LDV2T,M-Q) 233*> If JOBV2T = 'Y', V2T contains the (M-Q)-by-(M-Q) unitary 234*> matrix V2**H. 235*> \endverbatim 236*> 237*> \param[in] LDV2T 238*> \verbatim 239*> LDV2T is INTEGER 240*> The leading dimension of V2T. If JOBV2T = 'Y', LDV2T >= 241*> MAX(1,M-Q). 242*> \endverbatim 243*> 244*> \param[out] WORK 245*> \verbatim 246*> WORK is COMPLEX array, dimension (MAX(1,LWORK)) 247*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 248*> \endverbatim 249*> 250*> \param[in] LWORK 251*> \verbatim 252*> LWORK is INTEGER 253*> The dimension of the array WORK. 254*> 255*> If LWORK = -1, then a workspace query is assumed; the routine 256*> only calculates the optimal size of the WORK array, returns 257*> this value as the first entry of the work array, and no error 258*> message related to LWORK is issued by XERBLA. 259*> \endverbatim 260*> 261*> \param[out] RWORK 262*> \verbatim 263*> RWORK is REAL array, dimension MAX(1,LRWORK) 264*> On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK. 265*> If INFO > 0 on exit, RWORK(2:R) contains the values PHI(1), 266*> ..., PHI(R-1) that, together with THETA(1), ..., THETA(R), 267*> define the matrix in intermediate bidiagonal-block form 268*> remaining after nonconvergence. INFO specifies the number 269*> of nonzero PHI's. 270*> \endverbatim 271*> 272*> \param[in] LRWORK 273*> \verbatim 274*> LRWORK is INTEGER 275*> The dimension of the array RWORK. 276*> 277*> If LRWORK = -1, then a workspace query is assumed; the routine 278*> only calculates the optimal size of the RWORK array, returns 279*> this value as the first entry of the work array, and no error 280*> message related to LRWORK is issued by XERBLA. 281*> \endverbatim 282*> 283*> \param[out] IWORK 284*> \verbatim 285*> IWORK is INTEGER array, dimension (M-MIN(P,M-P,Q,M-Q)) 286*> \endverbatim 287*> 288*> \param[out] INFO 289*> \verbatim 290*> INFO is INTEGER 291*> = 0: successful exit. 292*> < 0: if INFO = -i, the i-th argument had an illegal value. 293*> > 0: CBBCSD did not converge. See the description of RWORK 294*> above for details. 295*> \endverbatim 296* 297*> \par References: 298* ================ 299*> 300*> [1] Brian D. Sutton. Computing the complete CS decomposition. Numer. 301*> Algorithms, 50(1):33-65, 2009. 302* 303* Authors: 304* ======== 305* 306*> \author Univ. of Tennessee 307*> \author Univ. of California Berkeley 308*> \author Univ. of Colorado Denver 309*> \author NAG Ltd. 310* 311*> \ingroup complexOTHERcomputational 312* 313* ===================================================================== 314 RECURSIVE SUBROUTINE CUNCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, 315 $ SIGNS, M, P, Q, X11, LDX11, X12, 316 $ LDX12, X21, LDX21, X22, LDX22, THETA, 317 $ U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, 318 $ LDV2T, WORK, LWORK, RWORK, LRWORK, 319 $ IWORK, INFO ) 320* 321* -- LAPACK computational routine -- 322* -- LAPACK is a software package provided by Univ. of Tennessee, -- 323* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 324* 325* .. Scalar Arguments .. 326 CHARACTER JOBU1, JOBU2, JOBV1T, JOBV2T, SIGNS, TRANS 327 INTEGER INFO, LDU1, LDU2, LDV1T, LDV2T, LDX11, LDX12, 328 $ LDX21, LDX22, LRWORK, LWORK, M, P, Q 329* .. 330* .. Array Arguments .. 331 INTEGER IWORK( * ) 332 REAL THETA( * ) 333 REAL RWORK( * ) 334 COMPLEX U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ), 335 $ V2T( LDV2T, * ), WORK( * ), X11( LDX11, * ), 336 $ X12( LDX12, * ), X21( LDX21, * ), X22( LDX22, 337 $ * ) 338* .. 339* 340* =================================================================== 341* 342* .. Parameters .. 343 COMPLEX ONE, ZERO 344 PARAMETER ( ONE = (1.0E0,0.0E0), 345 $ ZERO = (0.0E0,0.0E0) ) 346* .. 347* .. Local Scalars .. 348 CHARACTER TRANST, SIGNST 349 INTEGER CHILDINFO, I, IB11D, IB11E, IB12D, IB12E, 350 $ IB21D, IB21E, IB22D, IB22E, IBBCSD, IORBDB, 351 $ IORGLQ, IORGQR, IPHI, ITAUP1, ITAUP2, ITAUQ1, 352 $ ITAUQ2, J, LBBCSDWORK, LBBCSDWORKMIN, 353 $ LBBCSDWORKOPT, LORBDBWORK, LORBDBWORKMIN, 354 $ LORBDBWORKOPT, LORGLQWORK, LORGLQWORKMIN, 355 $ LORGLQWORKOPT, LORGQRWORK, LORGQRWORKMIN, 356 $ LORGQRWORKOPT, LWORKMIN, LWORKOPT, P1, Q1 357 LOGICAL COLMAJOR, DEFAULTSIGNS, LQUERY, WANTU1, WANTU2, 358 $ WANTV1T, WANTV2T 359 INTEGER LRWORKMIN, LRWORKOPT 360 LOGICAL LRQUERY 361* .. 362* .. External Subroutines .. 363 EXTERNAL XERBLA, CBBCSD, CLACPY, CLAPMR, CLAPMT, 364 $ CUNBDB, CUNGLQ, CUNGQR 365* .. 366* .. External Functions .. 367 LOGICAL LSAME 368 EXTERNAL LSAME 369* .. 370* .. Intrinsic Functions 371 INTRINSIC INT, MAX, MIN 372* .. 373* .. Executable Statements .. 374* 375* Test input arguments 376* 377 INFO = 0 378 WANTU1 = LSAME( JOBU1, 'Y' ) 379 WANTU2 = LSAME( JOBU2, 'Y' ) 380 WANTV1T = LSAME( JOBV1T, 'Y' ) 381 WANTV2T = LSAME( JOBV2T, 'Y' ) 382 COLMAJOR = .NOT. LSAME( TRANS, 'T' ) 383 DEFAULTSIGNS = .NOT. LSAME( SIGNS, 'O' ) 384 LQUERY = LWORK .EQ. -1 385 LRQUERY = LRWORK .EQ. -1 386 IF( M .LT. 0 ) THEN 387 INFO = -7 388 ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN 389 INFO = -8 390 ELSE IF( Q .LT. 0 .OR. Q .GT. M ) THEN 391 INFO = -9 392 ELSE IF ( COLMAJOR .AND. LDX11 .LT. MAX( 1, P ) ) THEN 393 INFO = -11 394 ELSE IF (.NOT. COLMAJOR .AND. LDX11 .LT. MAX( 1, Q ) ) THEN 395 INFO = -11 396 ELSE IF (COLMAJOR .AND. LDX12 .LT. MAX( 1, P ) ) THEN 397 INFO = -13 398 ELSE IF (.NOT. COLMAJOR .AND. LDX12 .LT. MAX( 1, M-Q ) ) THEN 399 INFO = -13 400 ELSE IF (COLMAJOR .AND. LDX21 .LT. MAX( 1, M-P ) ) THEN 401 INFO = -15 402 ELSE IF (.NOT. COLMAJOR .AND. LDX21 .LT. MAX( 1, Q ) ) THEN 403 INFO = -15 404 ELSE IF (COLMAJOR .AND. LDX22 .LT. MAX( 1, M-P ) ) THEN 405 INFO = -17 406 ELSE IF (.NOT. COLMAJOR .AND. LDX22 .LT. MAX( 1, M-Q ) ) THEN 407 INFO = -17 408 ELSE IF( WANTU1 .AND. LDU1 .LT. P ) THEN 409 INFO = -20 410 ELSE IF( WANTU2 .AND. LDU2 .LT. M-P ) THEN 411 INFO = -22 412 ELSE IF( WANTV1T .AND. LDV1T .LT. Q ) THEN 413 INFO = -24 414 ELSE IF( WANTV2T .AND. LDV2T .LT. M-Q ) THEN 415 INFO = -26 416 END IF 417* 418* Work with transpose if convenient 419* 420 IF( INFO .EQ. 0 .AND. MIN( P, M-P ) .LT. MIN( Q, M-Q ) ) THEN 421 IF( COLMAJOR ) THEN 422 TRANST = 'T' 423 ELSE 424 TRANST = 'N' 425 END IF 426 IF( DEFAULTSIGNS ) THEN 427 SIGNST = 'O' 428 ELSE 429 SIGNST = 'D' 430 END IF 431 CALL CUNCSD( JOBV1T, JOBV2T, JOBU1, JOBU2, TRANST, SIGNST, M, 432 $ Q, P, X11, LDX11, X21, LDX21, X12, LDX12, X22, 433 $ LDX22, THETA, V1T, LDV1T, V2T, LDV2T, U1, LDU1, 434 $ U2, LDU2, WORK, LWORK, RWORK, LRWORK, IWORK, 435 $ INFO ) 436 RETURN 437 END IF 438* 439* Work with permutation [ 0 I; I 0 ] * X * [ 0 I; I 0 ] if 440* convenient 441* 442 IF( INFO .EQ. 0 .AND. M-Q .LT. Q ) THEN 443 IF( DEFAULTSIGNS ) THEN 444 SIGNST = 'O' 445 ELSE 446 SIGNST = 'D' 447 END IF 448 CALL CUNCSD( JOBU2, JOBU1, JOBV2T, JOBV1T, TRANS, SIGNST, M, 449 $ M-P, M-Q, X22, LDX22, X21, LDX21, X12, LDX12, X11, 450 $ LDX11, THETA, U2, LDU2, U1, LDU1, V2T, LDV2T, V1T, 451 $ LDV1T, WORK, LWORK, RWORK, LRWORK, IWORK, INFO ) 452 RETURN 453 END IF 454* 455* Compute workspace 456* 457 IF( INFO .EQ. 0 ) THEN 458* 459* Real workspace 460* 461 IPHI = 2 462 IB11D = IPHI + MAX( 1, Q - 1 ) 463 IB11E = IB11D + MAX( 1, Q ) 464 IB12D = IB11E + MAX( 1, Q - 1 ) 465 IB12E = IB12D + MAX( 1, Q ) 466 IB21D = IB12E + MAX( 1, Q - 1 ) 467 IB21E = IB21D + MAX( 1, Q ) 468 IB22D = IB21E + MAX( 1, Q - 1 ) 469 IB22E = IB22D + MAX( 1, Q ) 470 IBBCSD = IB22E + MAX( 1, Q - 1 ) 471 CALL CBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, 472 $ THETA, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, 473 $ V2T, LDV2T, THETA, THETA, THETA, THETA, THETA, 474 $ THETA, THETA, THETA, RWORK, -1, CHILDINFO ) 475 LBBCSDWORKOPT = INT( RWORK(1) ) 476 LBBCSDWORKMIN = LBBCSDWORKOPT 477 LRWORKOPT = IBBCSD + LBBCSDWORKOPT - 1 478 LRWORKMIN = IBBCSD + LBBCSDWORKMIN - 1 479 RWORK(1) = LRWORKOPT 480* 481* Complex workspace 482* 483 ITAUP1 = 2 484 ITAUP2 = ITAUP1 + MAX( 1, P ) 485 ITAUQ1 = ITAUP2 + MAX( 1, M - P ) 486 ITAUQ2 = ITAUQ1 + MAX( 1, Q ) 487 IORGQR = ITAUQ2 + MAX( 1, M - Q ) 488 CALL CUNGQR( M-Q, M-Q, M-Q, U1, MAX(1,M-Q), U1, WORK, -1, 489 $ CHILDINFO ) 490 LORGQRWORKOPT = INT( WORK(1) ) 491 LORGQRWORKMIN = MAX( 1, M - Q ) 492 IORGLQ = ITAUQ2 + MAX( 1, M - Q ) 493 CALL CUNGLQ( M-Q, M-Q, M-Q, U1, MAX(1,M-Q), U1, WORK, -1, 494 $ CHILDINFO ) 495 LORGLQWORKOPT = INT( WORK(1) ) 496 LORGLQWORKMIN = MAX( 1, M - Q ) 497 IORBDB = ITAUQ2 + MAX( 1, M - Q ) 498 CALL CUNBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, 499 $ X21, LDX21, X22, LDX22, THETA, THETA, U1, U2, 500 $ V1T, V2T, WORK, -1, CHILDINFO ) 501 LORBDBWORKOPT = INT( WORK(1) ) 502 LORBDBWORKMIN = LORBDBWORKOPT 503 LWORKOPT = MAX( IORGQR + LORGQRWORKOPT, IORGLQ + LORGLQWORKOPT, 504 $ IORBDB + LORBDBWORKOPT ) - 1 505 LWORKMIN = MAX( IORGQR + LORGQRWORKMIN, IORGLQ + LORGLQWORKMIN, 506 $ IORBDB + LORBDBWORKMIN ) - 1 507 WORK(1) = MAX(LWORKOPT,LWORKMIN) 508* 509 IF( LWORK .LT. LWORKMIN 510 $ .AND. .NOT. ( LQUERY .OR. LRQUERY ) ) THEN 511 INFO = -22 512 ELSE IF( LRWORK .LT. LRWORKMIN 513 $ .AND. .NOT. ( LQUERY .OR. LRQUERY ) ) THEN 514 INFO = -24 515 ELSE 516 LORGQRWORK = LWORK - IORGQR + 1 517 LORGLQWORK = LWORK - IORGLQ + 1 518 LORBDBWORK = LWORK - IORBDB + 1 519 LBBCSDWORK = LRWORK - IBBCSD + 1 520 END IF 521 END IF 522* 523* Abort if any illegal arguments 524* 525 IF( INFO .NE. 0 ) THEN 526 CALL XERBLA( 'CUNCSD', -INFO ) 527 RETURN 528 ELSE IF( LQUERY .OR. LRQUERY ) THEN 529 RETURN 530 END IF 531* 532* Transform to bidiagonal block form 533* 534 CALL CUNBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, X21, 535 $ LDX21, X22, LDX22, THETA, RWORK(IPHI), WORK(ITAUP1), 536 $ WORK(ITAUP2), WORK(ITAUQ1), WORK(ITAUQ2), 537 $ WORK(IORBDB), LORBDBWORK, CHILDINFO ) 538* 539* Accumulate Householder reflectors 540* 541 IF( COLMAJOR ) THEN 542 IF( WANTU1 .AND. P .GT. 0 ) THEN 543 CALL CLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 ) 544 CALL CUNGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR), 545 $ LORGQRWORK, INFO) 546 END IF 547 IF( WANTU2 .AND. M-P .GT. 0 ) THEN 548 CALL CLACPY( 'L', M-P, Q, X21, LDX21, U2, LDU2 ) 549 CALL CUNGQR( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2), 550 $ WORK(IORGQR), LORGQRWORK, INFO ) 551 END IF 552 IF( WANTV1T .AND. Q .GT. 0 ) THEN 553 CALL CLACPY( 'U', Q-1, Q-1, X11(1,2), LDX11, V1T(2,2), 554 $ LDV1T ) 555 V1T(1, 1) = ONE 556 DO J = 2, Q 557 V1T(1,J) = ZERO 558 V1T(J,1) = ZERO 559 END DO 560 CALL CUNGLQ( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1), 561 $ WORK(IORGLQ), LORGLQWORK, INFO ) 562 END IF 563 IF( WANTV2T .AND. M-Q .GT. 0 ) THEN 564 CALL CLACPY( 'U', P, M-Q, X12, LDX12, V2T, LDV2T ) 565 IF( M-P .GT. Q ) THEN 566 CALL CLACPY( 'U', M-P-Q, M-P-Q, X22(Q+1,P+1), LDX22, 567 $ V2T(P+1,P+1), LDV2T ) 568 END IF 569 IF( M .GT. Q ) THEN 570 CALL CUNGLQ( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2), 571 $ WORK(IORGLQ), LORGLQWORK, INFO ) 572 END IF 573 END IF 574 ELSE 575 IF( WANTU1 .AND. P .GT. 0 ) THEN 576 CALL CLACPY( 'U', Q, P, X11, LDX11, U1, LDU1 ) 577 CALL CUNGLQ( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGLQ), 578 $ LORGLQWORK, INFO) 579 END IF 580 IF( WANTU2 .AND. M-P .GT. 0 ) THEN 581 CALL CLACPY( 'U', Q, M-P, X21, LDX21, U2, LDU2 ) 582 CALL CUNGLQ( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2), 583 $ WORK(IORGLQ), LORGLQWORK, INFO ) 584 END IF 585 IF( WANTV1T .AND. Q .GT. 0 ) THEN 586 CALL CLACPY( 'L', Q-1, Q-1, X11(2,1), LDX11, V1T(2,2), 587 $ LDV1T ) 588 V1T(1, 1) = ONE 589 DO J = 2, Q 590 V1T(1,J) = ZERO 591 V1T(J,1) = ZERO 592 END DO 593 CALL CUNGQR( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1), 594 $ WORK(IORGQR), LORGQRWORK, INFO ) 595 END IF 596 IF( WANTV2T .AND. M-Q .GT. 0 ) THEN 597 P1 = MIN( P+1, M ) 598 Q1 = MIN( Q+1, M ) 599 CALL CLACPY( 'L', M-Q, P, X12, LDX12, V2T, LDV2T ) 600 IF ( M .GT. P+Q ) THEN 601 CALL CLACPY( 'L', M-P-Q, M-P-Q, X22(P1,Q1), LDX22, 602 $ V2T(P+1,P+1), LDV2T ) 603 END IF 604 CALL CUNGQR( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2), 605 $ WORK(IORGQR), LORGQRWORK, INFO ) 606 END IF 607 END IF 608* 609* Compute the CSD of the matrix in bidiagonal-block form 610* 611 CALL CBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, THETA, 612 $ RWORK(IPHI), U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, 613 $ LDV2T, RWORK(IB11D), RWORK(IB11E), RWORK(IB12D), 614 $ RWORK(IB12E), RWORK(IB21D), RWORK(IB21E), 615 $ RWORK(IB22D), RWORK(IB22E), RWORK(IBBCSD), 616 $ LBBCSDWORK, INFO ) 617* 618* Permute rows and columns to place identity submatrices in top- 619* left corner of (1,1)-block and/or bottom-right corner of (1,2)- 620* block and/or bottom-right corner of (2,1)-block and/or top-left 621* corner of (2,2)-block 622* 623 IF( Q .GT. 0 .AND. WANTU2 ) THEN 624 DO I = 1, Q 625 IWORK(I) = M - P - Q + I 626 END DO 627 DO I = Q + 1, M - P 628 IWORK(I) = I - Q 629 END DO 630 IF( COLMAJOR ) THEN 631 CALL CLAPMT( .FALSE., M-P, M-P, U2, LDU2, IWORK ) 632 ELSE 633 CALL CLAPMR( .FALSE., M-P, M-P, U2, LDU2, IWORK ) 634 END IF 635 END IF 636 IF( M .GT. 0 .AND. WANTV2T ) THEN 637 DO I = 1, P 638 IWORK(I) = M - P - Q + I 639 END DO 640 DO I = P + 1, M - Q 641 IWORK(I) = I - P 642 END DO 643 IF( .NOT. COLMAJOR ) THEN 644 CALL CLAPMT( .FALSE., M-Q, M-Q, V2T, LDV2T, IWORK ) 645 ELSE 646 CALL CLAPMR( .FALSE., M-Q, M-Q, V2T, LDV2T, IWORK ) 647 END IF 648 END IF 649* 650 RETURN 651* 652* End CUNCSD 653* 654 END 655 656