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