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