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*> \date June 2016 312* 313*> \ingroup complexOTHERcomputational 314* 315* ===================================================================== 316 RECURSIVE SUBROUTINE CUNCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, 317 $ SIGNS, M, P, Q, X11, LDX11, X12, 318 $ LDX12, X21, LDX21, X22, LDX22, THETA, 319 $ U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, 320 $ LDV2T, WORK, LWORK, RWORK, LRWORK, 321 $ IWORK, INFO ) 322* 323* -- LAPACK computational routine (version 3.7.1) -- 324* -- LAPACK is a software package provided by Univ. of Tennessee, -- 325* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 326* June 2016 327* 328* .. Scalar Arguments .. 329 CHARACTER JOBU1, JOBU2, JOBV1T, JOBV2T, SIGNS, TRANS 330 INTEGER INFO, LDU1, LDU2, LDV1T, LDV2T, LDX11, LDX12, 331 $ LDX21, LDX22, LRWORK, LWORK, M, P, Q 332* .. 333* .. Array Arguments .. 334 INTEGER IWORK( * ) 335 REAL THETA( * ) 336 REAL RWORK( * ) 337 COMPLEX U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ), 338 $ V2T( LDV2T, * ), WORK( * ), X11( LDX11, * ), 339 $ X12( LDX12, * ), X21( LDX21, * ), X22( LDX22, 340 $ * ) 341* .. 342* 343* =================================================================== 344* 345* .. Parameters .. 346 COMPLEX ONE, ZERO 347 PARAMETER ( ONE = (1.0E0,0.0E0), 348 $ ZERO = (0.0E0,0.0E0) ) 349* .. 350* .. Local Scalars .. 351 CHARACTER TRANST, SIGNST 352 INTEGER CHILDINFO, I, IB11D, IB11E, IB12D, IB12E, 353 $ IB21D, IB21E, IB22D, IB22E, IBBCSD, IORBDB, 354 $ IORGLQ, IORGQR, IPHI, ITAUP1, ITAUP2, ITAUQ1, 355 $ ITAUQ2, J, LBBCSDWORK, LBBCSDWORKMIN, 356 $ LBBCSDWORKOPT, LORBDBWORK, LORBDBWORKMIN, 357 $ LORBDBWORKOPT, LORGLQWORK, LORGLQWORKMIN, 358 $ LORGLQWORKOPT, LORGQRWORK, LORGQRWORKMIN, 359 $ LORGQRWORKOPT, LWORKMIN, LWORKOPT, P1, Q1 360 LOGICAL COLMAJOR, DEFAULTSIGNS, LQUERY, WANTU1, WANTU2, 361 $ WANTV1T, WANTV2T 362 INTEGER LRWORKMIN, LRWORKOPT 363 LOGICAL LRQUERY 364* .. 365* .. External Subroutines .. 366 EXTERNAL XERBLA, CBBCSD, CLACPY, CLAPMR, CLAPMT, 367 $ CUNBDB, CUNGLQ, CUNGQR 368* .. 369* .. External Functions .. 370 LOGICAL LSAME 371 EXTERNAL LSAME 372* .. 373* .. Intrinsic Functions 374 INTRINSIC INT, MAX, MIN 375* .. 376* .. Executable Statements .. 377* 378* Test input arguments 379* 380 INFO = 0 381 WANTU1 = LSAME( JOBU1, 'Y' ) 382 WANTU2 = LSAME( JOBU2, 'Y' ) 383 WANTV1T = LSAME( JOBV1T, 'Y' ) 384 WANTV2T = LSAME( JOBV2T, 'Y' ) 385 COLMAJOR = .NOT. LSAME( TRANS, 'T' ) 386 DEFAULTSIGNS = .NOT. LSAME( SIGNS, 'O' ) 387 LQUERY = LWORK .EQ. -1 388 LRQUERY = LRWORK .EQ. -1 389 IF( M .LT. 0 ) THEN 390 INFO = -7 391 ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN 392 INFO = -8 393 ELSE IF( Q .LT. 0 .OR. Q .GT. M ) THEN 394 INFO = -9 395 ELSE IF ( COLMAJOR .AND. LDX11 .LT. MAX( 1, P ) ) THEN 396 INFO = -11 397 ELSE IF (.NOT. COLMAJOR .AND. LDX11 .LT. MAX( 1, Q ) ) THEN 398 INFO = -11 399 ELSE IF (COLMAJOR .AND. LDX12 .LT. MAX( 1, P ) ) THEN 400 INFO = -13 401 ELSE IF (.NOT. COLMAJOR .AND. LDX12 .LT. MAX( 1, M-Q ) ) THEN 402 INFO = -13 403 ELSE IF (COLMAJOR .AND. LDX21 .LT. MAX( 1, M-P ) ) THEN 404 INFO = -15 405 ELSE IF (.NOT. COLMAJOR .AND. LDX21 .LT. MAX( 1, Q ) ) THEN 406 INFO = -15 407 ELSE IF (COLMAJOR .AND. LDX22 .LT. MAX( 1, M-P ) ) THEN 408 INFO = -17 409 ELSE IF (.NOT. COLMAJOR .AND. LDX22 .LT. MAX( 1, M-Q ) ) THEN 410 INFO = -17 411 ELSE IF( WANTU1 .AND. LDU1 .LT. P ) THEN 412 INFO = -20 413 ELSE IF( WANTU2 .AND. LDU2 .LT. M-P ) THEN 414 INFO = -22 415 ELSE IF( WANTV1T .AND. LDV1T .LT. Q ) THEN 416 INFO = -24 417 ELSE IF( WANTV2T .AND. LDV2T .LT. M-Q ) THEN 418 INFO = -26 419 END IF 420* 421* Work with transpose if convenient 422* 423 IF( INFO .EQ. 0 .AND. MIN( P, M-P ) .LT. MIN( Q, M-Q ) ) THEN 424 IF( COLMAJOR ) THEN 425 TRANST = 'T' 426 ELSE 427 TRANST = 'N' 428 END IF 429 IF( DEFAULTSIGNS ) THEN 430 SIGNST = 'O' 431 ELSE 432 SIGNST = 'D' 433 END IF 434 CALL CUNCSD( JOBV1T, JOBV2T, JOBU1, JOBU2, TRANST, SIGNST, M, 435 $ Q, P, X11, LDX11, X21, LDX21, X12, LDX12, X22, 436 $ LDX22, THETA, V1T, LDV1T, V2T, LDV2T, U1, LDU1, 437 $ U2, LDU2, WORK, LWORK, RWORK, LRWORK, IWORK, 438 $ INFO ) 439 RETURN 440 END IF 441* 442* Work with permutation [ 0 I; I 0 ] * X * [ 0 I; I 0 ] if 443* convenient 444* 445 IF( INFO .EQ. 0 .AND. M-Q .LT. Q ) THEN 446 IF( DEFAULTSIGNS ) THEN 447 SIGNST = 'O' 448 ELSE 449 SIGNST = 'D' 450 END IF 451 CALL CUNCSD( JOBU2, JOBU1, JOBV2T, JOBV1T, TRANS, SIGNST, M, 452 $ M-P, M-Q, X22, LDX22, X21, LDX21, X12, LDX12, X11, 453 $ LDX11, THETA, U2, LDU2, U1, LDU1, V2T, LDV2T, V1T, 454 $ LDV1T, WORK, LWORK, RWORK, LRWORK, IWORK, INFO ) 455 RETURN 456 END IF 457* 458* Compute workspace 459* 460 IF( INFO .EQ. 0 ) THEN 461* 462* Real workspace 463* 464 IPHI = 2 465 IB11D = IPHI + MAX( 1, Q - 1 ) 466 IB11E = IB11D + MAX( 1, Q ) 467 IB12D = IB11E + MAX( 1, Q - 1 ) 468 IB12E = IB12D + MAX( 1, Q ) 469 IB21D = IB12E + MAX( 1, Q - 1 ) 470 IB21E = IB21D + MAX( 1, Q ) 471 IB22D = IB21E + MAX( 1, Q - 1 ) 472 IB22E = IB22D + MAX( 1, Q ) 473 IBBCSD = IB22E + MAX( 1, Q - 1 ) 474 CALL CBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, 475 $ THETA, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, 476 $ V2T, LDV2T, THETA, THETA, THETA, THETA, THETA, 477 $ THETA, THETA, THETA, RWORK, -1, CHILDINFO ) 478 LBBCSDWORKOPT = INT( RWORK(1) ) 479 LBBCSDWORKMIN = LBBCSDWORKOPT 480 LRWORKOPT = IBBCSD + LBBCSDWORKOPT - 1 481 LRWORKMIN = IBBCSD + LBBCSDWORKMIN - 1 482 RWORK(1) = LRWORKOPT 483* 484* Complex workspace 485* 486 ITAUP1 = 2 487 ITAUP2 = ITAUP1 + MAX( 1, P ) 488 ITAUQ1 = ITAUP2 + MAX( 1, M - P ) 489 ITAUQ2 = ITAUQ1 + MAX( 1, Q ) 490 IORGQR = ITAUQ2 + MAX( 1, M - Q ) 491 CALL CUNGQR( M-Q, M-Q, M-Q, U1, MAX(1,M-Q), U1, WORK, -1, 492 $ CHILDINFO ) 493 LORGQRWORKOPT = INT( WORK(1) ) 494 LORGQRWORKMIN = MAX( 1, M - Q ) 495 IORGLQ = ITAUQ2 + MAX( 1, M - Q ) 496 CALL CUNGLQ( M-Q, M-Q, M-Q, U1, MAX(1,M-Q), U1, WORK, -1, 497 $ CHILDINFO ) 498 LORGLQWORKOPT = INT( WORK(1) ) 499 LORGLQWORKMIN = MAX( 1, M - Q ) 500 IORBDB = ITAUQ2 + MAX( 1, M - Q ) 501 CALL CUNBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, 502 $ X21, LDX21, X22, LDX22, THETA, THETA, U1, U2, 503 $ V1T, V2T, WORK, -1, CHILDINFO ) 504 LORBDBWORKOPT = INT( WORK(1) ) 505 LORBDBWORKMIN = LORBDBWORKOPT 506 LWORKOPT = MAX( IORGQR + LORGQRWORKOPT, IORGLQ + LORGLQWORKOPT, 507 $ IORBDB + LORBDBWORKOPT ) - 1 508 LWORKMIN = MAX( IORGQR + LORGQRWORKMIN, IORGLQ + LORGLQWORKMIN, 509 $ IORBDB + LORBDBWORKMIN ) - 1 510 WORK(1) = MAX(LWORKOPT,LWORKMIN) 511* 512 IF( LWORK .LT. LWORKMIN 513 $ .AND. .NOT. ( LQUERY .OR. LRQUERY ) ) THEN 514 INFO = -22 515 ELSE IF( LRWORK .LT. LRWORKMIN 516 $ .AND. .NOT. ( LQUERY .OR. LRQUERY ) ) THEN 517 INFO = -24 518 ELSE 519 LORGQRWORK = LWORK - IORGQR + 1 520 LORGLQWORK = LWORK - IORGLQ + 1 521 LORBDBWORK = LWORK - IORBDB + 1 522 LBBCSDWORK = LRWORK - IBBCSD + 1 523 END IF 524 END IF 525* 526* Abort if any illegal arguments 527* 528 IF( INFO .NE. 0 ) THEN 529 CALL XERBLA( 'CUNCSD', -INFO ) 530 RETURN 531 ELSE IF( LQUERY .OR. LRQUERY ) THEN 532 RETURN 533 END IF 534* 535* Transform to bidiagonal block form 536* 537 CALL CUNBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, X21, 538 $ LDX21, X22, LDX22, THETA, RWORK(IPHI), WORK(ITAUP1), 539 $ WORK(ITAUP2), WORK(ITAUQ1), WORK(ITAUQ2), 540 $ WORK(IORBDB), LORBDBWORK, CHILDINFO ) 541* 542* Accumulate Householder reflectors 543* 544 IF( COLMAJOR ) THEN 545 IF( WANTU1 .AND. P .GT. 0 ) THEN 546 CALL CLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 ) 547 CALL CUNGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR), 548 $ LORGQRWORK, INFO) 549 END IF 550 IF( WANTU2 .AND. M-P .GT. 0 ) THEN 551 CALL CLACPY( 'L', M-P, Q, X21, LDX21, U2, LDU2 ) 552 CALL CUNGQR( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2), 553 $ WORK(IORGQR), LORGQRWORK, INFO ) 554 END IF 555 IF( WANTV1T .AND. Q .GT. 0 ) THEN 556 CALL CLACPY( 'U', Q-1, Q-1, X11(1,2), LDX11, V1T(2,2), 557 $ LDV1T ) 558 V1T(1, 1) = ONE 559 DO J = 2, Q 560 V1T(1,J) = ZERO 561 V1T(J,1) = ZERO 562 END DO 563 CALL CUNGLQ( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1), 564 $ WORK(IORGLQ), LORGLQWORK, INFO ) 565 END IF 566 IF( WANTV2T .AND. M-Q .GT. 0 ) THEN 567 CALL CLACPY( 'U', P, M-Q, X12, LDX12, V2T, LDV2T ) 568 IF( M-P .GT. Q ) THEN 569 CALL CLACPY( 'U', M-P-Q, M-P-Q, X22(Q+1,P+1), LDX22, 570 $ V2T(P+1,P+1), LDV2T ) 571 END IF 572 IF( M .GT. Q ) THEN 573 CALL CUNGLQ( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2), 574 $ WORK(IORGLQ), LORGLQWORK, INFO ) 575 END IF 576 END IF 577 ELSE 578 IF( WANTU1 .AND. P .GT. 0 ) THEN 579 CALL CLACPY( 'U', Q, P, X11, LDX11, U1, LDU1 ) 580 CALL CUNGLQ( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGLQ), 581 $ LORGLQWORK, INFO) 582 END IF 583 IF( WANTU2 .AND. M-P .GT. 0 ) THEN 584 CALL CLACPY( 'U', Q, M-P, X21, LDX21, U2, LDU2 ) 585 CALL CUNGLQ( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2), 586 $ WORK(IORGLQ), LORGLQWORK, INFO ) 587 END IF 588 IF( WANTV1T .AND. Q .GT. 0 ) THEN 589 CALL CLACPY( 'L', Q-1, Q-1, X11(2,1), LDX11, V1T(2,2), 590 $ LDV1T ) 591 V1T(1, 1) = ONE 592 DO J = 2, Q 593 V1T(1,J) = ZERO 594 V1T(J,1) = ZERO 595 END DO 596 CALL CUNGQR( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1), 597 $ WORK(IORGQR), LORGQRWORK, INFO ) 598 END IF 599 IF( WANTV2T .AND. M-Q .GT. 0 ) THEN 600 P1 = MIN( P+1, M ) 601 Q1 = MIN( Q+1, M ) 602 CALL CLACPY( 'L', M-Q, P, X12, LDX12, V2T, LDV2T ) 603 IF ( M .GT. P+Q ) THEN 604 CALL CLACPY( 'L', M-P-Q, M-P-Q, X22(P1,Q1), LDX22, 605 $ V2T(P+1,P+1), LDV2T ) 606 END IF 607 CALL CUNGQR( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2), 608 $ WORK(IORGQR), LORGQRWORK, INFO ) 609 END IF 610 END IF 611* 612* Compute the CSD of the matrix in bidiagonal-block form 613* 614 CALL CBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, THETA, 615 $ RWORK(IPHI), U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, 616 $ LDV2T, RWORK(IB11D), RWORK(IB11E), RWORK(IB12D), 617 $ RWORK(IB12E), RWORK(IB21D), RWORK(IB21E), 618 $ RWORK(IB22D), RWORK(IB22E), RWORK(IBBCSD), 619 $ LBBCSDWORK, INFO ) 620* 621* Permute rows and columns to place identity submatrices in top- 622* left corner of (1,1)-block and/or bottom-right corner of (1,2)- 623* block and/or bottom-right corner of (2,1)-block and/or top-left 624* corner of (2,2)-block 625* 626 IF( Q .GT. 0 .AND. WANTU2 ) THEN 627 DO I = 1, Q 628 IWORK(I) = M - P - Q + I 629 END DO 630 DO I = Q + 1, M - P 631 IWORK(I) = I - Q 632 END DO 633 IF( COLMAJOR ) THEN 634 CALL CLAPMT( .FALSE., M-P, M-P, U2, LDU2, IWORK ) 635 ELSE 636 CALL CLAPMR( .FALSE., M-P, M-P, U2, LDU2, IWORK ) 637 END IF 638 END IF 639 IF( M .GT. 0 .AND. WANTV2T ) THEN 640 DO I = 1, P 641 IWORK(I) = M - P - Q + I 642 END DO 643 DO I = P + 1, M - Q 644 IWORK(I) = I - P 645 END DO 646 IF( .NOT. COLMAJOR ) THEN 647 CALL CLAPMT( .FALSE., M-Q, M-Q, V2T, LDV2T, IWORK ) 648 ELSE 649 CALL CLAPMR( .FALSE., M-Q, M-Q, V2T, LDV2T, IWORK ) 650 END IF 651 END IF 652* 653 RETURN 654* 655* End CUNCSD 656* 657 END 658 659