1*> \brief \b ZUNCSD2BY1 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download ZUNCSD2BY1 + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zuncsd2by1.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zuncsd2by1.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zuncsd2by1.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE ZUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, 22* X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T, 23* LDV1T, WORK, LWORK, RWORK, LRWORK, IWORK, 24* INFO ) 25* 26* .. Scalar Arguments .. 27* CHARACTER JOBU1, JOBU2, JOBV1T 28* INTEGER INFO, LDU1, LDU2, LDV1T, LWORK, LDX11, LDX21, 29* $ M, P, Q 30* INTEGER LRWORK, LRWORKMIN, LRWORKOPT 31* .. 32* .. Array Arguments .. 33* DOUBLE PRECISION RWORK(*) 34* DOUBLE PRECISION THETA(*) 35* COMPLEX*16 U1(LDU1,*), U2(LDU2,*), V1T(LDV1T,*), WORK(*), 36* $ X11(LDX11,*), X21(LDX21,*) 37* INTEGER IWORK(*) 38* .. 39* 40* 41*> \par Purpose: 42* ============= 43*> 44*>\verbatim 45*> 46*> ZUNCSD2BY1 computes the CS decomposition of an M-by-Q matrix X with 47*> orthonormal columns that has been partitioned into a 2-by-1 block 48*> structure: 49*> 50*> [ I1 0 0 ] 51*> [ 0 C 0 ] 52*> [ X11 ] [ U1 | ] [ 0 0 0 ] 53*> X = [-----] = [---------] [----------] V1**T . 54*> [ X21 ] [ | U2 ] [ 0 0 0 ] 55*> [ 0 S 0 ] 56*> [ 0 0 I2] 57*> 58*> X11 is P-by-Q. The unitary matrices U1, U2, and V1 are P-by-P, 59*> (M-P)-by-(M-P), and Q-by-Q, respectively. C and S are R-by-R 60*> nonnegative diagonal matrices satisfying C^2 + S^2 = I, in which 61*> R = MIN(P,M-P,Q,M-Q). I1 is a K1-by-K1 identity matrix and I2 is a 62*> K2-by-K2 identity matrix, where K1 = MAX(Q+P-M,0), K2 = MAX(Q-P,0). 63*> \endverbatim 64* 65* Arguments: 66* ========== 67* 68*> \param[in] JOBU1 69*> \verbatim 70*> JOBU1 is CHARACTER 71*> = 'Y': U1 is computed; 72*> otherwise: U1 is not computed. 73*> \endverbatim 74*> 75*> \param[in] JOBU2 76*> \verbatim 77*> JOBU2 is CHARACTER 78*> = 'Y': U2 is computed; 79*> otherwise: U2 is not computed. 80*> \endverbatim 81*> 82*> \param[in] JOBV1T 83*> \verbatim 84*> JOBV1T is CHARACTER 85*> = 'Y': V1T is computed; 86*> otherwise: V1T is not computed. 87*> \endverbatim 88*> 89*> \param[in] M 90*> \verbatim 91*> M is INTEGER 92*> The number of rows in X. 93*> \endverbatim 94*> 95*> \param[in] P 96*> \verbatim 97*> P is INTEGER 98*> The number of rows in X11. 0 <= P <= M. 99*> \endverbatim 100*> 101*> \param[in] Q 102*> \verbatim 103*> Q is INTEGER 104*> The number of columns in X11 and X21. 0 <= Q <= M. 105*> \endverbatim 106*> 107*> \param[in,out] X11 108*> \verbatim 109*> X11 is COMPLEX*16 array, dimension (LDX11,Q) 110*> On entry, part of the unitary matrix whose CSD is desired. 111*> \endverbatim 112*> 113*> \param[in] LDX11 114*> \verbatim 115*> LDX11 is INTEGER 116*> The leading dimension of X11. LDX11 >= MAX(1,P). 117*> \endverbatim 118*> 119*> \param[in,out] X21 120*> \verbatim 121*> X21 is COMPLEX*16 array, dimension (LDX21,Q) 122*> On entry, part of the unitary matrix whose CSD is desired. 123*> \endverbatim 124*> 125*> \param[in] LDX21 126*> \verbatim 127*> LDX21 is INTEGER 128*> The leading dimension of X21. LDX21 >= MAX(1,M-P). 129*> \endverbatim 130*> 131*> \param[out] THETA 132*> \verbatim 133*> THETA is DOUBLE PRECISION array, dimension (R), in which R = 134*> MIN(P,M-P,Q,M-Q). 135*> C = DIAG( COS(THETA(1)), ... , COS(THETA(R)) ) and 136*> S = DIAG( SIN(THETA(1)), ... , SIN(THETA(R)) ). 137*> \endverbatim 138*> 139*> \param[out] U1 140*> \verbatim 141*> U1 is COMPLEX*16 array, dimension (P) 142*> If JOBU1 = 'Y', U1 contains the P-by-P unitary matrix U1. 143*> \endverbatim 144*> 145*> \param[in] LDU1 146*> \verbatim 147*> LDU1 is INTEGER 148*> The leading dimension of U1. If JOBU1 = 'Y', LDU1 >= 149*> MAX(1,P). 150*> \endverbatim 151*> 152*> \param[out] U2 153*> \verbatim 154*> U2 is COMPLEX*16 array, dimension (M-P) 155*> If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) unitary 156*> matrix U2. 157*> \endverbatim 158*> 159*> \param[in] LDU2 160*> \verbatim 161*> LDU2 is INTEGER 162*> The leading dimension of U2. If JOBU2 = 'Y', LDU2 >= 163*> MAX(1,M-P). 164*> \endverbatim 165*> 166*> \param[out] V1T 167*> \verbatim 168*> V1T is COMPLEX*16 array, dimension (Q) 169*> If JOBV1T = 'Y', V1T contains the Q-by-Q matrix unitary 170*> matrix V1**T. 171*> \endverbatim 172*> 173*> \param[in] LDV1T 174*> \verbatim 175*> LDV1T is INTEGER 176*> The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >= 177*> MAX(1,Q). 178*> \endverbatim 179*> 180*> \param[out] WORK 181*> \verbatim 182*> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK)) 183*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 184*> \endverbatim 185*> 186*> \param[in] LWORK 187*> \verbatim 188*> LWORK is INTEGER 189*> The dimension of the array WORK. 190*> 191*> If LWORK = -1, then a workspace query is assumed; the routine 192*> only calculates the optimal size of the WORK and RWORK 193*> arrays, returns this value as the first entry of the WORK 194*> and RWORK array, respectively, and no error message related 195*> to LWORK or LRWORK is issued by XERBLA. 196*> \endverbatim 197*> 198*> \param[out] RWORK 199*> \verbatim 200*> RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK)) 201*> On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK. 202*> If INFO > 0 on exit, RWORK(2:R) contains the values PHI(1), 203*> ..., PHI(R-1) that, together with THETA(1), ..., THETA(R), 204*> define the matrix in intermediate bidiagonal-block form 205*> remaining after nonconvergence. INFO specifies the number 206*> of nonzero PHI's. 207*> \endverbatim 208*> 209*> \param[in] LRWORK 210*> \verbatim 211*> LRWORK is INTEGER 212*> The dimension of the array RWORK. 213*> 214*> If LRWORK=-1, then a workspace query is assumed; the routine 215*> only calculates the optimal size of the WORK and RWORK 216*> arrays, returns this value as the first entry of the WORK 217*> and RWORK array, respectively, and no error message related 218*> to LWORK or LRWORK is issued by XERBLA. 219*> \endverbatim 220* 221*> \param[out] IWORK 222*> \verbatim 223*> IWORK is INTEGER array, dimension (M-MIN(P,M-P,Q,M-Q)) 224*> \endverbatim 225*> 226*> \param[out] INFO 227*> \verbatim 228*> INFO is INTEGER 229*> = 0: successful exit. 230*> < 0: if INFO = -i, the i-th argument had an illegal value. 231*> > 0: ZBBCSD did not converge. See the description of WORK 232*> above for details. 233*> \endverbatim 234* 235*> \par References: 236* ================ 237*> 238*> [1] Brian D. Sutton. Computing the complete CS decomposition. Numer. 239*> Algorithms, 50(1):33-65, 2009. 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*> \ingroup complex16OTHERcomputational 250* 251* ===================================================================== 252 SUBROUTINE ZUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, 253 $ X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T, 254 $ LDV1T, WORK, LWORK, RWORK, LRWORK, IWORK, 255 $ INFO ) 256* 257* -- LAPACK computational routine -- 258* -- LAPACK is a software package provided by Univ. of Tennessee, -- 259* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 260* 261* .. Scalar Arguments .. 262 CHARACTER JOBU1, JOBU2, JOBV1T 263 INTEGER INFO, LDU1, LDU2, LDV1T, LWORK, LDX11, LDX21, 264 $ M, P, Q 265 INTEGER LRWORK, LRWORKMIN, LRWORKOPT 266* .. 267* .. Array Arguments .. 268 DOUBLE PRECISION RWORK(*) 269 DOUBLE PRECISION THETA(*) 270 COMPLEX*16 U1(LDU1,*), U2(LDU2,*), V1T(LDV1T,*), WORK(*), 271 $ X11(LDX11,*), X21(LDX21,*) 272 INTEGER IWORK(*) 273* .. 274* 275* ===================================================================== 276* 277* .. Parameters .. 278 COMPLEX*16 ONE, ZERO 279 PARAMETER ( ONE = (1.0D0,0.0D0), ZERO = (0.0D0,0.0D0) ) 280* .. 281* .. Local Scalars .. 282 INTEGER CHILDINFO, I, IB11D, IB11E, IB12D, IB12E, 283 $ IB21D, IB21E, IB22D, IB22E, IBBCSD, IORBDB, 284 $ IORGLQ, IORGQR, IPHI, ITAUP1, ITAUP2, ITAUQ1, 285 $ J, LBBCSD, LORBDB, LORGLQ, LORGLQMIN, 286 $ LORGLQOPT, LORGQR, LORGQRMIN, LORGQROPT, 287 $ LWORKMIN, LWORKOPT, R 288 LOGICAL LQUERY, WANTU1, WANTU2, WANTV1T 289* .. 290* .. Local Arrays .. 291 DOUBLE PRECISION DUM( 1 ) 292 COMPLEX*16 CDUM( 1, 1 ) 293* .. 294* .. External Subroutines .. 295 EXTERNAL ZBBCSD, ZCOPY, ZLACPY, ZLAPMR, ZLAPMT, ZUNBDB1, 296 $ ZUNBDB2, ZUNBDB3, ZUNBDB4, ZUNGLQ, ZUNGQR, 297 $ XERBLA 298* .. 299* .. External Functions .. 300 LOGICAL LSAME 301 EXTERNAL LSAME 302* .. 303* .. Intrinsic Function .. 304 INTRINSIC INT, MAX, MIN 305* .. 306* .. Executable Statements .. 307* 308* Test input arguments 309* 310 INFO = 0 311 WANTU1 = LSAME( JOBU1, 'Y' ) 312 WANTU2 = LSAME( JOBU2, 'Y' ) 313 WANTV1T = LSAME( JOBV1T, 'Y' ) 314 LQUERY = ( LWORK.EQ.-1 ) .OR. ( LRWORK.EQ.-1 ) 315* 316 IF( M .LT. 0 ) THEN 317 INFO = -4 318 ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN 319 INFO = -5 320 ELSE IF( Q .LT. 0 .OR. Q .GT. M ) THEN 321 INFO = -6 322 ELSE IF( LDX11 .LT. MAX( 1, P ) ) THEN 323 INFO = -8 324 ELSE IF( LDX21 .LT. MAX( 1, M-P ) ) THEN 325 INFO = -10 326 ELSE IF( WANTU1 .AND. LDU1 .LT. MAX( 1, P ) ) THEN 327 INFO = -13 328 ELSE IF( WANTU2 .AND. LDU2 .LT. MAX( 1, M - P ) ) THEN 329 INFO = -15 330 ELSE IF( WANTV1T .AND. LDV1T .LT. MAX( 1, Q ) ) THEN 331 INFO = -17 332 END IF 333* 334 R = MIN( P, M-P, Q, M-Q ) 335* 336* Compute workspace 337* 338* WORK layout: 339* |-----------------------------------------| 340* | LWORKOPT (1) | 341* |-----------------------------------------| 342* | TAUP1 (MAX(1,P)) | 343* | TAUP2 (MAX(1,M-P)) | 344* | TAUQ1 (MAX(1,Q)) | 345* |-----------------------------------------| 346* | ZUNBDB WORK | ZUNGQR WORK | ZUNGLQ WORK | 347* | | | | 348* | | | | 349* | | | | 350* | | | | 351* |-----------------------------------------| 352* RWORK layout: 353* |------------------| 354* | LRWORKOPT (1) | 355* |------------------| 356* | PHI (MAX(1,R-1)) | 357* |------------------| 358* | B11D (R) | 359* | B11E (R-1) | 360* | B12D (R) | 361* | B12E (R-1) | 362* | B21D (R) | 363* | B21E (R-1) | 364* | B22D (R) | 365* | B22E (R-1) | 366* | ZBBCSD RWORK | 367* |------------------| 368* 369 IF( INFO .EQ. 0 ) THEN 370 IPHI = 2 371 IB11D = IPHI + MAX( 1, R-1 ) 372 IB11E = IB11D + MAX( 1, R ) 373 IB12D = IB11E + MAX( 1, R - 1 ) 374 IB12E = IB12D + MAX( 1, R ) 375 IB21D = IB12E + MAX( 1, R - 1 ) 376 IB21E = IB21D + MAX( 1, R ) 377 IB22D = IB21E + MAX( 1, R - 1 ) 378 IB22E = IB22D + MAX( 1, R ) 379 IBBCSD = IB22E + MAX( 1, R - 1 ) 380 ITAUP1 = 2 381 ITAUP2 = ITAUP1 + MAX( 1, P ) 382 ITAUQ1 = ITAUP2 + MAX( 1, M-P ) 383 IORBDB = ITAUQ1 + MAX( 1, Q ) 384 IORGQR = ITAUQ1 + MAX( 1, Q ) 385 IORGLQ = ITAUQ1 + MAX( 1, Q ) 386 LORGQRMIN = 1 387 LORGQROPT = 1 388 LORGLQMIN = 1 389 LORGLQOPT = 1 390 IF( R .EQ. Q ) THEN 391 CALL ZUNBDB1( M, P, Q, X11, LDX11, X21, LDX21, THETA, DUM, 392 $ CDUM, CDUM, CDUM, WORK, -1, CHILDINFO ) 393 LORBDB = INT( WORK(1) ) 394 IF( WANTU1 .AND. P .GT. 0 ) THEN 395 CALL ZUNGQR( P, P, Q, U1, LDU1, CDUM, WORK(1), -1, 396 $ CHILDINFO ) 397 LORGQRMIN = MAX( LORGQRMIN, P ) 398 LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) ) 399 ENDIF 400 IF( WANTU2 .AND. M-P .GT. 0 ) THEN 401 CALL ZUNGQR( M-P, M-P, Q, U2, LDU2, CDUM, WORK(1), -1, 402 $ CHILDINFO ) 403 LORGQRMIN = MAX( LORGQRMIN, M-P ) 404 LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) ) 405 END IF 406 IF( WANTV1T .AND. Q .GT. 0 ) THEN 407 CALL ZUNGLQ( Q-1, Q-1, Q-1, V1T, LDV1T, 408 $ CDUM, WORK(1), -1, CHILDINFO ) 409 LORGLQMIN = MAX( LORGLQMIN, Q-1 ) 410 LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) ) 411 END IF 412 CALL ZBBCSD( JOBU1, JOBU2, JOBV1T, 'N', 'N', M, P, Q, THETA, 413 $ DUM, U1, LDU1, U2, LDU2, V1T, LDV1T, CDUM, 1, 414 $ DUM, DUM, DUM, DUM, DUM, DUM, DUM, DUM, 415 $ RWORK(1), -1, CHILDINFO ) 416 LBBCSD = INT( RWORK(1) ) 417 ELSE IF( R .EQ. P ) THEN 418 CALL ZUNBDB2( M, P, Q, X11, LDX11, X21, LDX21, THETA, DUM, 419 $ CDUM, CDUM, CDUM, WORK(1), -1, CHILDINFO ) 420 LORBDB = INT( WORK(1) ) 421 IF( WANTU1 .AND. P .GT. 0 ) THEN 422 CALL ZUNGQR( P-1, P-1, P-1, U1(2,2), LDU1, CDUM, WORK(1), 423 $ -1, CHILDINFO ) 424 LORGQRMIN = MAX( LORGQRMIN, P-1 ) 425 LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) ) 426 END IF 427 IF( WANTU2 .AND. M-P .GT. 0 ) THEN 428 CALL ZUNGQR( M-P, M-P, Q, U2, LDU2, CDUM, WORK(1), -1, 429 $ CHILDINFO ) 430 LORGQRMIN = MAX( LORGQRMIN, M-P ) 431 LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) ) 432 END IF 433 IF( WANTV1T .AND. Q .GT. 0 ) THEN 434 CALL ZUNGLQ( Q, Q, R, V1T, LDV1T, CDUM, WORK(1), -1, 435 $ CHILDINFO ) 436 LORGLQMIN = MAX( LORGLQMIN, Q ) 437 LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) ) 438 END IF 439 CALL ZBBCSD( JOBV1T, 'N', JOBU1, JOBU2, 'T', M, Q, P, THETA, 440 $ DUM, V1T, LDV1T, CDUM, 1, U1, LDU1, U2, LDU2, 441 $ DUM, DUM, DUM, DUM, DUM, DUM, DUM, DUM, 442 $ RWORK(1), -1, CHILDINFO ) 443 LBBCSD = INT( RWORK(1) ) 444 ELSE IF( R .EQ. M-P ) THEN 445 CALL ZUNBDB3( M, P, Q, X11, LDX11, X21, LDX21, THETA, DUM, 446 $ CDUM, CDUM, CDUM, WORK(1), -1, CHILDINFO ) 447 LORBDB = INT( WORK(1) ) 448 IF( WANTU1 .AND. P .GT. 0 ) THEN 449 CALL ZUNGQR( P, P, Q, U1, LDU1, CDUM, WORK(1), -1, 450 $ CHILDINFO ) 451 LORGQRMIN = MAX( LORGQRMIN, P ) 452 LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) ) 453 END IF 454 IF( WANTU2 .AND. M-P .GT. 0 ) THEN 455 CALL ZUNGQR( M-P-1, M-P-1, M-P-1, U2(2,2), LDU2, CDUM, 456 $ WORK(1), -1, CHILDINFO ) 457 LORGQRMIN = MAX( LORGQRMIN, M-P-1 ) 458 LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) ) 459 END IF 460 IF( WANTV1T .AND. Q .GT. 0 ) THEN 461 CALL ZUNGLQ( Q, Q, R, V1T, LDV1T, CDUM, WORK(1), -1, 462 $ CHILDINFO ) 463 LORGLQMIN = MAX( LORGLQMIN, Q ) 464 LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) ) 465 END IF 466 CALL ZBBCSD( 'N', JOBV1T, JOBU2, JOBU1, 'T', M, M-Q, M-P, 467 $ THETA, DUM, CDUM, 1, V1T, LDV1T, U2, LDU2, U1, 468 $ LDU1, DUM, DUM, DUM, DUM, DUM, DUM, DUM, DUM, 469 $ RWORK(1), -1, CHILDINFO ) 470 LBBCSD = INT( RWORK(1) ) 471 ELSE 472 CALL ZUNBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA, DUM, 473 $ CDUM, CDUM, CDUM, CDUM, WORK(1), -1, CHILDINFO 474 $ ) 475 LORBDB = M + INT( WORK(1) ) 476 IF( WANTU1 .AND. P .GT. 0 ) THEN 477 CALL ZUNGQR( P, P, M-Q, U1, LDU1, CDUM, WORK(1), -1, 478 $ CHILDINFO ) 479 LORGQRMIN = MAX( LORGQRMIN, P ) 480 LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) ) 481 END IF 482 IF( WANTU2 .AND. M-P .GT. 0 ) THEN 483 CALL ZUNGQR( M-P, M-P, M-Q, U2, LDU2, CDUM, WORK(1), -1, 484 $ CHILDINFO ) 485 LORGQRMIN = MAX( LORGQRMIN, M-P ) 486 LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) ) 487 END IF 488 IF( WANTV1T .AND. Q .GT. 0 ) THEN 489 CALL ZUNGLQ( Q, Q, Q, V1T, LDV1T, CDUM, WORK(1), -1, 490 $ CHILDINFO ) 491 LORGLQMIN = MAX( LORGLQMIN, Q ) 492 LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) ) 493 END IF 494 CALL ZBBCSD( JOBU2, JOBU1, 'N', JOBV1T, 'N', M, M-P, M-Q, 495 $ THETA, DUM, U2, LDU2, U1, LDU1, CDUM, 1, V1T, 496 $ LDV1T, DUM, DUM, DUM, DUM, DUM, DUM, DUM, DUM, 497 $ RWORK(1), -1, CHILDINFO ) 498 LBBCSD = INT( RWORK(1) ) 499 END IF 500 LRWORKMIN = IBBCSD+LBBCSD-1 501 LRWORKOPT = LRWORKMIN 502 RWORK(1) = LRWORKOPT 503 LWORKMIN = MAX( IORBDB+LORBDB-1, 504 $ IORGQR+LORGQRMIN-1, 505 $ IORGLQ+LORGLQMIN-1 ) 506 LWORKOPT = MAX( IORBDB+LORBDB-1, 507 $ IORGQR+LORGQROPT-1, 508 $ IORGLQ+LORGLQOPT-1 ) 509 WORK(1) = LWORKOPT 510 IF( LWORK .LT. LWORKMIN .AND. .NOT.LQUERY ) THEN 511 INFO = -19 512 END IF 513 IF( LRWORK .LT. LRWORKMIN .AND. .NOT.LQUERY ) THEN 514 INFO = -21 515 END IF 516 END IF 517 IF( INFO .NE. 0 ) THEN 518 CALL XERBLA( 'ZUNCSD2BY1', -INFO ) 519 RETURN 520 ELSE IF( LQUERY ) THEN 521 RETURN 522 END IF 523 LORGQR = LWORK-IORGQR+1 524 LORGLQ = LWORK-IORGLQ+1 525* 526* Handle four cases separately: R = Q, R = P, R = M-P, and R = M-Q, 527* in which R = MIN(P,M-P,Q,M-Q) 528* 529 IF( R .EQ. Q ) THEN 530* 531* Case 1: R = Q 532* 533* Simultaneously bidiagonalize X11 and X21 534* 535 CALL ZUNBDB1( M, P, Q, X11, LDX11, X21, LDX21, THETA, 536 $ RWORK(IPHI), WORK(ITAUP1), WORK(ITAUP2), 537 $ WORK(ITAUQ1), WORK(IORBDB), LORBDB, CHILDINFO ) 538* 539* Accumulate Householder reflectors 540* 541 IF( WANTU1 .AND. P .GT. 0 ) THEN 542 CALL ZLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 ) 543 CALL ZUNGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR), 544 $ LORGQR, CHILDINFO ) 545 END IF 546 IF( WANTU2 .AND. M-P .GT. 0 ) THEN 547 CALL ZLACPY( 'L', M-P, Q, X21, LDX21, U2, LDU2 ) 548 CALL ZUNGQR( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2), 549 $ WORK(IORGQR), LORGQR, CHILDINFO ) 550 END IF 551 IF( WANTV1T .AND. Q .GT. 0 ) THEN 552 V1T(1,1) = ONE 553 DO J = 2, Q 554 V1T(1,J) = ZERO 555 V1T(J,1) = ZERO 556 END DO 557 CALL ZLACPY( 'U', Q-1, Q-1, X21(1,2), LDX21, V1T(2,2), 558 $ LDV1T ) 559 CALL ZUNGLQ( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1), 560 $ WORK(IORGLQ), LORGLQ, CHILDINFO ) 561 END IF 562* 563* Simultaneously diagonalize X11 and X21. 564* 565 CALL ZBBCSD( JOBU1, JOBU2, JOBV1T, 'N', 'N', M, P, Q, THETA, 566 $ RWORK(IPHI), U1, LDU1, U2, LDU2, V1T, LDV1T, CDUM, 567 $ 1, RWORK(IB11D), RWORK(IB11E), RWORK(IB12D), 568 $ RWORK(IB12E), RWORK(IB21D), RWORK(IB21E), 569 $ RWORK(IB22D), RWORK(IB22E), RWORK(IBBCSD), 570 $ LRWORK-IBBCSD+1, CHILDINFO ) 571* 572* Permute rows and columns to place zero submatrices in 573* preferred positions 574* 575 IF( Q .GT. 0 .AND. WANTU2 ) THEN 576 DO I = 1, Q 577 IWORK(I) = M - P - Q + I 578 END DO 579 DO I = Q + 1, M - P 580 IWORK(I) = I - Q 581 END DO 582 CALL ZLAPMT( .FALSE., M-P, M-P, U2, LDU2, IWORK ) 583 END IF 584 ELSE IF( R .EQ. P ) THEN 585* 586* Case 2: R = P 587* 588* Simultaneously bidiagonalize X11 and X21 589* 590 CALL ZUNBDB2( M, P, Q, X11, LDX11, X21, LDX21, THETA, 591 $ RWORK(IPHI), WORK(ITAUP1), WORK(ITAUP2), 592 $ WORK(ITAUQ1), WORK(IORBDB), LORBDB, CHILDINFO ) 593* 594* Accumulate Householder reflectors 595* 596 IF( WANTU1 .AND. P .GT. 0 ) THEN 597 U1(1,1) = ONE 598 DO J = 2, P 599 U1(1,J) = ZERO 600 U1(J,1) = ZERO 601 END DO 602 CALL ZLACPY( 'L', P-1, P-1, X11(2,1), LDX11, U1(2,2), LDU1 ) 603 CALL ZUNGQR( P-1, P-1, P-1, U1(2,2), LDU1, WORK(ITAUP1), 604 $ WORK(IORGQR), LORGQR, CHILDINFO ) 605 END IF 606 IF( WANTU2 .AND. M-P .GT. 0 ) THEN 607 CALL ZLACPY( 'L', M-P, Q, X21, LDX21, U2, LDU2 ) 608 CALL ZUNGQR( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2), 609 $ WORK(IORGQR), LORGQR, CHILDINFO ) 610 END IF 611 IF( WANTV1T .AND. Q .GT. 0 ) THEN 612 CALL ZLACPY( 'U', P, Q, X11, LDX11, V1T, LDV1T ) 613 CALL ZUNGLQ( Q, Q, R, V1T, LDV1T, WORK(ITAUQ1), 614 $ WORK(IORGLQ), LORGLQ, CHILDINFO ) 615 END IF 616* 617* Simultaneously diagonalize X11 and X21. 618* 619 CALL ZBBCSD( JOBV1T, 'N', JOBU1, JOBU2, 'T', M, Q, P, THETA, 620 $ RWORK(IPHI), V1T, LDV1T, CDUM, 1, U1, LDU1, U2, 621 $ LDU2, RWORK(IB11D), RWORK(IB11E), RWORK(IB12D), 622 $ RWORK(IB12E), RWORK(IB21D), RWORK(IB21E), 623 $ RWORK(IB22D), RWORK(IB22E), RWORK(IBBCSD), LBBCSD, 624 $ CHILDINFO ) 625* 626* Permute rows and columns to place identity submatrices in 627* preferred positions 628* 629 IF( Q .GT. 0 .AND. WANTU2 ) THEN 630 DO I = 1, Q 631 IWORK(I) = M - P - Q + I 632 END DO 633 DO I = Q + 1, M - P 634 IWORK(I) = I - Q 635 END DO 636 CALL ZLAPMT( .FALSE., M-P, M-P, U2, LDU2, IWORK ) 637 END IF 638 ELSE IF( R .EQ. M-P ) THEN 639* 640* Case 3: R = M-P 641* 642* Simultaneously bidiagonalize X11 and X21 643* 644 CALL ZUNBDB3( M, P, Q, X11, LDX11, X21, LDX21, THETA, 645 $ RWORK(IPHI), WORK(ITAUP1), WORK(ITAUP2), 646 $ WORK(ITAUQ1), WORK(IORBDB), LORBDB, CHILDINFO ) 647* 648* Accumulate Householder reflectors 649* 650 IF( WANTU1 .AND. P .GT. 0 ) THEN 651 CALL ZLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 ) 652 CALL ZUNGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR), 653 $ LORGQR, CHILDINFO ) 654 END IF 655 IF( WANTU2 .AND. M-P .GT. 0 ) THEN 656 U2(1,1) = ONE 657 DO J = 2, M-P 658 U2(1,J) = ZERO 659 U2(J,1) = ZERO 660 END DO 661 CALL ZLACPY( 'L', M-P-1, M-P-1, X21(2,1), LDX21, U2(2,2), 662 $ LDU2 ) 663 CALL ZUNGQR( M-P-1, M-P-1, M-P-1, U2(2,2), LDU2, 664 $ WORK(ITAUP2), WORK(IORGQR), LORGQR, CHILDINFO ) 665 END IF 666 IF( WANTV1T .AND. Q .GT. 0 ) THEN 667 CALL ZLACPY( 'U', M-P, Q, X21, LDX21, V1T, LDV1T ) 668 CALL ZUNGLQ( Q, Q, R, V1T, LDV1T, WORK(ITAUQ1), 669 $ WORK(IORGLQ), LORGLQ, CHILDINFO ) 670 END IF 671* 672* Simultaneously diagonalize X11 and X21. 673* 674 CALL ZBBCSD( 'N', JOBV1T, JOBU2, JOBU1, 'T', M, M-Q, M-P, 675 $ THETA, RWORK(IPHI), CDUM, 1, V1T, LDV1T, U2, LDU2, 676 $ U1, LDU1, RWORK(IB11D), RWORK(IB11E), 677 $ RWORK(IB12D), RWORK(IB12E), RWORK(IB21D), 678 $ RWORK(IB21E), RWORK(IB22D), RWORK(IB22E), 679 $ RWORK(IBBCSD), LBBCSD, CHILDINFO ) 680* 681* Permute rows and columns to place identity submatrices in 682* preferred positions 683* 684 IF( Q .GT. R ) THEN 685 DO I = 1, R 686 IWORK(I) = Q - R + I 687 END DO 688 DO I = R + 1, Q 689 IWORK(I) = I - R 690 END DO 691 IF( WANTU1 ) THEN 692 CALL ZLAPMT( .FALSE., P, Q, U1, LDU1, IWORK ) 693 END IF 694 IF( WANTV1T ) THEN 695 CALL ZLAPMR( .FALSE., Q, Q, V1T, LDV1T, IWORK ) 696 END IF 697 END IF 698 ELSE 699* 700* Case 4: R = M-Q 701* 702* Simultaneously bidiagonalize X11 and X21 703* 704 CALL ZUNBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA, 705 $ RWORK(IPHI), WORK(ITAUP1), WORK(ITAUP2), 706 $ WORK(ITAUQ1), WORK(IORBDB), WORK(IORBDB+M), 707 $ LORBDB-M, CHILDINFO ) 708* 709* Accumulate Householder reflectors 710* 711 IF( WANTU2 .AND. M-P .GT. 0 ) THEN 712 CALL ZCOPY( M-P, WORK(IORBDB+P), 1, U2, 1 ) 713 END IF 714 IF( WANTU1 .AND. P .GT. 0 ) THEN 715 CALL ZCOPY( P, WORK(IORBDB), 1, U1, 1 ) 716 DO J = 2, P 717 U1(1,J) = ZERO 718 END DO 719 CALL ZLACPY( 'L', P-1, M-Q-1, X11(2,1), LDX11, U1(2,2), 720 $ LDU1 ) 721 CALL ZUNGQR( P, P, M-Q, U1, LDU1, WORK(ITAUP1), 722 $ WORK(IORGQR), LORGQR, CHILDINFO ) 723 END IF 724 IF( WANTU2 .AND. M-P .GT. 0 ) THEN 725 DO J = 2, M-P 726 U2(1,J) = ZERO 727 END DO 728 CALL ZLACPY( 'L', M-P-1, M-Q-1, X21(2,1), LDX21, U2(2,2), 729 $ LDU2 ) 730 CALL ZUNGQR( M-P, M-P, M-Q, U2, LDU2, WORK(ITAUP2), 731 $ WORK(IORGQR), LORGQR, CHILDINFO ) 732 END IF 733 IF( WANTV1T .AND. Q .GT. 0 ) THEN 734 CALL ZLACPY( 'U', M-Q, Q, X21, LDX21, V1T, LDV1T ) 735 CALL ZLACPY( 'U', P-(M-Q), Q-(M-Q), X11(M-Q+1,M-Q+1), LDX11, 736 $ V1T(M-Q+1,M-Q+1), LDV1T ) 737 CALL ZLACPY( 'U', -P+Q, Q-P, X21(M-Q+1,P+1), LDX21, 738 $ V1T(P+1,P+1), LDV1T ) 739 CALL ZUNGLQ( Q, Q, Q, V1T, LDV1T, WORK(ITAUQ1), 740 $ WORK(IORGLQ), LORGLQ, CHILDINFO ) 741 END IF 742* 743* Simultaneously diagonalize X11 and X21. 744* 745 CALL ZBBCSD( JOBU2, JOBU1, 'N', JOBV1T, 'N', M, M-P, M-Q, 746 $ THETA, RWORK(IPHI), U2, LDU2, U1, LDU1, CDUM, 1, 747 $ V1T, LDV1T, RWORK(IB11D), RWORK(IB11E), 748 $ RWORK(IB12D), RWORK(IB12E), RWORK(IB21D), 749 $ RWORK(IB21E), RWORK(IB22D), RWORK(IB22E), 750 $ RWORK(IBBCSD), LBBCSD, CHILDINFO ) 751* 752* Permute rows and columns to place identity submatrices in 753* preferred positions 754* 755 IF( P .GT. R ) THEN 756 DO I = 1, R 757 IWORK(I) = P - R + I 758 END DO 759 DO I = R + 1, P 760 IWORK(I) = I - R 761 END DO 762 IF( WANTU1 ) THEN 763 CALL ZLAPMT( .FALSE., P, P, U1, LDU1, IWORK ) 764 END IF 765 IF( WANTV1T ) THEN 766 CALL ZLAPMR( .FALSE., P, Q, V1T, LDV1T, IWORK ) 767 END IF 768 END IF 769 END IF 770* 771 RETURN 772* 773* End of ZUNCSD2BY1 774* 775 END 776 777