![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DORBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, 2: $ X21, LDX21, X22, LDX22, THETA, PHI, TAUP1, 3: $ TAUP2, TAUQ1, TAUQ2, WORK, LWORK, INFO ) 4: IMPLICIT NONE 5: * 6: * -- LAPACK routine (version 3.3.0) -- 7: * 8: * -- Contributed by Brian Sutton of the Randolph-Macon College -- 9: * -- November 2010 10: * 11: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 12: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 13: * 14: * .. Scalar Arguments .. 15: CHARACTER SIGNS, TRANS 16: INTEGER INFO, LDX11, LDX12, LDX21, LDX22, LWORK, M, P, 17: $ Q 18: * .. 19: * .. Array Arguments .. 20: DOUBLE PRECISION PHI( * ), THETA( * ) 21: DOUBLE PRECISION TAUP1( * ), TAUP2( * ), TAUQ1( * ), TAUQ2( * ), 22: $ WORK( * ), X11( LDX11, * ), X12( LDX12, * ), 23: $ X21( LDX21, * ), X22( LDX22, * ) 24: * .. 25: * 26: * Purpose 27: * ======= 28: * 29: * DORBDB simultaneously bidiagonalizes the blocks of an M-by-M 30: * partitioned orthogonal matrix X: 31: * 32: * [ B11 | B12 0 0 ] 33: * [ X11 | X12 ] [ P1 | ] [ 0 | 0 -I 0 ] [ Q1 | ]**T 34: * X = [-----------] = [---------] [----------------] [---------] . 35: * [ X21 | X22 ] [ | P2 ] [ B21 | B22 0 0 ] [ | Q2 ] 36: * [ 0 | 0 0 I ] 37: * 38: * X11 is P-by-Q. Q must be no larger than P, M-P, or M-Q. (If this is 39: * not the case, then X must be transposed and/or permuted. This can be 40: * done in constant time using the TRANS and SIGNS options. See DORCSD 41: * for details.) 42: * 43: * The orthogonal matrices P1, P2, Q1, and Q2 are P-by-P, (M-P)-by- 44: * (M-P), Q-by-Q, and (M-Q)-by-(M-Q), respectively. They are 45: * represented implicitly by Householder vectors. 46: * 47: * B11, B12, B21, and B22 are Q-by-Q bidiagonal matrices represented 48: * implicitly by angles THETA, PHI. 49: * 50: * Arguments 51: * ========= 52: * 53: * TRANS (input) CHARACTER 54: * = 'T': X, U1, U2, V1T, and V2T are stored in row-major 55: * order; 56: * otherwise: X, U1, U2, V1T, and V2T are stored in column- 57: * major order. 58: * 59: * SIGNS (input) CHARACTER 60: * = 'O': The lower-left block is made nonpositive (the 61: * "other" convention); 62: * otherwise: The upper-right block is made nonpositive (the 63: * "default" convention). 64: * 65: * M (input) INTEGER 66: * The number of rows and columns in X. 67: * 68: * P (input) INTEGER 69: * The number of rows in X11 and X12. 0 <= P <= M. 70: * 71: * Q (input) INTEGER 72: * The number of columns in X11 and X21. 0 <= Q <= 73: * MIN(P,M-P,M-Q). 74: * 75: * X11 (input/output) DOUBLE PRECISION array, dimension (LDX11,Q) 76: * On entry, the top-left block of the orthogonal matrix to be 77: * reduced. On exit, the form depends on TRANS: 78: * If TRANS = 'N', then 79: * the columns of tril(X11) specify reflectors for P1, 80: * the rows of triu(X11,1) specify reflectors for Q1; 81: * else TRANS = 'T', and 82: * the rows of triu(X11) specify reflectors for P1, 83: * the columns of tril(X11,-1) specify reflectors for Q1. 84: * 85: * LDX11 (input) INTEGER 86: * The leading dimension of X11. If TRANS = 'N', then LDX11 >= 87: * P; else LDX11 >= Q. 88: * 89: * X12 (input/output) DOUBLE PRECISION array, dimension (LDX12,M-Q) 90: * On entry, the top-right block of the orthogonal matrix to 91: * be reduced. On exit, the form depends on TRANS: 92: * If TRANS = 'N', then 93: * the rows of triu(X12) specify the first P reflectors for 94: * Q2; 95: * else TRANS = 'T', and 96: * the columns of tril(X12) specify the first P reflectors 97: * for Q2. 98: * 99: * LDX12 (input) INTEGER 100: * The leading dimension of X12. If TRANS = 'N', then LDX12 >= 101: * P; else LDX11 >= M-Q. 102: * 103: * X21 (input/output) DOUBLE PRECISION array, dimension (LDX21,Q) 104: * On entry, the bottom-left block of the orthogonal matrix to 105: * be reduced. On exit, the form depends on TRANS: 106: * If TRANS = 'N', then 107: * the columns of tril(X21) specify reflectors for P2; 108: * else TRANS = 'T', and 109: * the rows of triu(X21) specify reflectors for P2. 110: * 111: * LDX21 (input) INTEGER 112: * The leading dimension of X21. If TRANS = 'N', then LDX21 >= 113: * M-P; else LDX21 >= Q. 114: * 115: * X22 (input/output) DOUBLE PRECISION array, dimension (LDX22,M-Q) 116: * On entry, the bottom-right block of the orthogonal matrix to 117: * be reduced. On exit, the form depends on TRANS: 118: * If TRANS = 'N', then 119: * the rows of triu(X22(Q+1:M-P,P+1:M-Q)) specify the last 120: * M-P-Q reflectors for Q2, 121: * else TRANS = 'T', and 122: * the columns of tril(X22(P+1:M-Q,Q+1:M-P)) specify the last 123: * M-P-Q reflectors for P2. 124: * 125: * LDX22 (input) INTEGER 126: * The leading dimension of X22. If TRANS = 'N', then LDX22 >= 127: * M-P; else LDX22 >= M-Q. 128: * 129: * THETA (output) DOUBLE PRECISION array, dimension (Q) 130: * The entries of the bidiagonal blocks B11, B12, B21, B22 can 131: * be computed from the angles THETA and PHI. See Further 132: * Details. 133: * 134: * PHI (output) DOUBLE PRECISION array, dimension (Q-1) 135: * The entries of the bidiagonal blocks B11, B12, B21, B22 can 136: * be computed from the angles THETA and PHI. See Further 137: * Details. 138: * 139: * TAUP1 (output) DOUBLE PRECISION array, dimension (P) 140: * The scalar factors of the elementary reflectors that define 141: * P1. 142: * 143: * TAUP2 (output) DOUBLE PRECISION array, dimension (M-P) 144: * The scalar factors of the elementary reflectors that define 145: * P2. 146: * 147: * TAUQ1 (output) DOUBLE PRECISION array, dimension (Q) 148: * The scalar factors of the elementary reflectors that define 149: * Q1. 150: * 151: * TAUQ2 (output) DOUBLE PRECISION array, dimension (M-Q) 152: * The scalar factors of the elementary reflectors that define 153: * Q2. 154: * 155: * WORK (workspace) DOUBLE PRECISION array, dimension (LWORK) 156: * 157: * LWORK (input) INTEGER 158: * The dimension of the array WORK. LWORK >= M-Q. 159: * 160: * If LWORK = -1, then a workspace query is assumed; the routine 161: * only calculates the optimal size of the WORK array, returns 162: * this value as the first entry of the WORK array, and no error 163: * message related to LWORK is issued by XERBLA. 164: * 165: * INFO (output) INTEGER 166: * = 0: successful exit. 167: * < 0: if INFO = -i, the i-th argument had an illegal value. 168: * 169: * Further Details 170: * =============== 171: * 172: * The bidiagonal blocks B11, B12, B21, and B22 are represented 173: * implicitly by angles THETA(1), ..., THETA(Q) and PHI(1), ..., 174: * PHI(Q-1). B11 and B21 are upper bidiagonal, while B21 and B22 are 175: * lower bidiagonal. Every entry in each bidiagonal band is a product 176: * of a sine or cosine of a THETA with a sine or cosine of a PHI. See 177: * [1] or DORCSD for details. 178: * 179: * P1, P2, Q1, and Q2 are represented as products of elementary 180: * reflectors. See DORCSD for details on generating P1, P2, Q1, and Q2 181: * using DORGQR and DORGLQ. 182: * 183: * Reference 184: * ========= 185: * 186: * [1] Brian D. Sutton. Computing the complete CS decomposition. Numer. 187: * Algorithms, 50(1):33-65, 2009. 188: * 189: * ==================================================================== 190: * 191: * .. Parameters .. 192: DOUBLE PRECISION REALONE 193: PARAMETER ( REALONE = 1.0D0 ) 194: DOUBLE PRECISION NEGONE, ONE 195: PARAMETER ( NEGONE = -1.0D0, ONE = 1.0D0 ) 196: * .. 197: * .. Local Scalars .. 198: LOGICAL COLMAJOR, LQUERY 199: INTEGER I, LWORKMIN, LWORKOPT 200: DOUBLE PRECISION Z1, Z2, Z3, Z4 201: * .. 202: * .. External Subroutines .. 203: EXTERNAL DAXPY, DLARF, DLARFGP, DSCAL, XERBLA 204: * .. 205: * .. External Functions .. 206: DOUBLE PRECISION DNRM2 207: LOGICAL LSAME 208: EXTERNAL DNRM2, LSAME 209: * .. 210: * .. Intrinsic Functions 211: INTRINSIC ATAN2, COS, MAX, MIN, SIN 212: * .. 213: * .. Executable Statements .. 214: * 215: * Test input arguments 216: * 217: INFO = 0 218: COLMAJOR = .NOT. LSAME( TRANS, 'T' ) 219: IF( .NOT. LSAME( SIGNS, 'O' ) ) THEN 220: Z1 = REALONE 221: Z2 = REALONE 222: Z3 = REALONE 223: Z4 = REALONE 224: ELSE 225: Z1 = REALONE 226: Z2 = -REALONE 227: Z3 = REALONE 228: Z4 = -REALONE 229: END IF 230: LQUERY = LWORK .EQ. -1 231: * 232: IF( M .LT. 0 ) THEN 233: INFO = -3 234: ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN 235: INFO = -4 236: ELSE IF( Q .LT. 0 .OR. Q .GT. P .OR. Q .GT. M-P .OR. 237: $ Q .GT. M-Q ) THEN 238: INFO = -5 239: ELSE IF( COLMAJOR .AND. LDX11 .LT. MAX( 1, P ) ) THEN 240: INFO = -7 241: ELSE IF( .NOT.COLMAJOR .AND. LDX11 .LT. MAX( 1, Q ) ) THEN 242: INFO = -7 243: ELSE IF( COLMAJOR .AND. LDX12 .LT. MAX( 1, P ) ) THEN 244: INFO = -9 245: ELSE IF( .NOT.COLMAJOR .AND. LDX12 .LT. MAX( 1, M-Q ) ) THEN 246: INFO = -9 247: ELSE IF( COLMAJOR .AND. LDX21 .LT. MAX( 1, M-P ) ) THEN 248: INFO = -11 249: ELSE IF( .NOT.COLMAJOR .AND. LDX21 .LT. MAX( 1, Q ) ) THEN 250: INFO = -11 251: ELSE IF( COLMAJOR .AND. LDX22 .LT. MAX( 1, M-P ) ) THEN 252: INFO = -13 253: ELSE IF( .NOT.COLMAJOR .AND. LDX22 .LT. MAX( 1, M-Q ) ) THEN 254: INFO = -13 255: END IF 256: * 257: * Compute workspace 258: * 259: IF( INFO .EQ. 0 ) THEN 260: LWORKOPT = M - Q 261: LWORKMIN = M - Q 262: WORK(1) = LWORKOPT 263: IF( LWORK .LT. LWORKMIN .AND. .NOT. LQUERY ) THEN 264: INFO = -21 265: END IF 266: END IF 267: IF( INFO .NE. 0 ) THEN 268: CALL XERBLA( 'xORBDB', -INFO ) 269: RETURN 270: ELSE IF( LQUERY ) THEN 271: RETURN 272: END IF 273: * 274: * Handle column-major and row-major separately 275: * 276: IF( COLMAJOR ) THEN 277: * 278: * Reduce columns 1, ..., Q of X11, X12, X21, and X22 279: * 280: DO I = 1, Q 281: * 282: IF( I .EQ. 1 ) THEN 283: CALL DSCAL( P-I+1, Z1, X11(I,I), 1 ) 284: ELSE 285: CALL DSCAL( P-I+1, Z1*COS(PHI(I-1)), X11(I,I), 1 ) 286: CALL DAXPY( P-I+1, -Z1*Z3*Z4*SIN(PHI(I-1)), X12(I,I-1), 287: $ 1, X11(I,I), 1 ) 288: END IF 289: IF( I .EQ. 1 ) THEN 290: CALL DSCAL( M-P-I+1, Z2, X21(I,I), 1 ) 291: ELSE 292: CALL DSCAL( M-P-I+1, Z2*COS(PHI(I-1)), X21(I,I), 1 ) 293: CALL DAXPY( M-P-I+1, -Z2*Z3*Z4*SIN(PHI(I-1)), X22(I,I-1), 294: $ 1, X21(I,I), 1 ) 295: END IF 296: * 297: THETA(I) = ATAN2( DNRM2( M-P-I+1, X21(I,I), 1 ), 298: $ DNRM2( P-I+1, X11(I,I), 1 ) ) 299: * 300: CALL DLARFGP( P-I+1, X11(I,I), X11(I+1,I), 1, TAUP1(I) ) 301: X11(I,I) = ONE 302: CALL DLARFGP( M-P-I+1, X21(I,I), X21(I+1,I), 1, TAUP2(I) ) 303: X21(I,I) = ONE 304: * 305: CALL DLARF( 'L', P-I+1, Q-I, X11(I,I), 1, TAUP1(I), 306: $ X11(I,I+1), LDX11, WORK ) 307: CALL DLARF( 'L', P-I+1, M-Q-I+1, X11(I,I), 1, TAUP1(I), 308: $ X12(I,I), LDX12, WORK ) 309: CALL DLARF( 'L', M-P-I+1, Q-I, X21(I,I), 1, TAUP2(I), 310: $ X21(I,I+1), LDX21, WORK ) 311: CALL DLARF( 'L', M-P-I+1, M-Q-I+1, X21(I,I), 1, TAUP2(I), 312: $ X22(I,I), LDX22, WORK ) 313: * 314: IF( I .LT. Q ) THEN 315: CALL DSCAL( Q-I, -Z1*Z3*SIN(THETA(I)), X11(I,I+1), 316: $ LDX11 ) 317: CALL DAXPY( Q-I, Z2*Z3*COS(THETA(I)), X21(I,I+1), LDX21, 318: $ X11(I,I+1), LDX11 ) 319: END IF 320: CALL DSCAL( M-Q-I+1, -Z1*Z4*SIN(THETA(I)), X12(I,I), LDX12 ) 321: CALL DAXPY( M-Q-I+1, Z2*Z4*COS(THETA(I)), X22(I,I), LDX22, 322: $ X12(I,I), LDX12 ) 323: * 324: IF( I .LT. Q ) 325: $ PHI(I) = ATAN2( DNRM2( Q-I, X11(I,I+1), LDX11 ), 326: $ DNRM2( M-Q-I+1, X12(I,I), LDX12 ) ) 327: * 328: IF( I .LT. Q ) THEN 329: CALL DLARFGP( Q-I, X11(I,I+1), X11(I,I+2), LDX11, 330: $ TAUQ1(I) ) 331: X11(I,I+1) = ONE 332: END IF 333: CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12, 334: $ TAUQ2(I) ) 335: X12(I,I) = ONE 336: * 337: IF( I .LT. Q ) THEN 338: CALL DLARF( 'R', P-I, Q-I, X11(I,I+1), LDX11, TAUQ1(I), 339: $ X11(I+1,I+1), LDX11, WORK ) 340: CALL DLARF( 'R', M-P-I, Q-I, X11(I,I+1), LDX11, TAUQ1(I), 341: $ X21(I+1,I+1), LDX21, WORK ) 342: END IF 343: CALL DLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I), 344: $ X12(I+1,I), LDX12, WORK ) 345: CALL DLARF( 'R', M-P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I), 346: $ X22(I+1,I), LDX22, WORK ) 347: * 348: END DO 349: * 350: * Reduce columns Q + 1, ..., P of X12, X22 351: * 352: DO I = Q + 1, P 353: * 354: CALL DSCAL( M-Q-I+1, -Z1*Z4, X12(I,I), LDX12 ) 355: CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12, 356: $ TAUQ2(I) ) 357: X12(I,I) = ONE 358: * 359: CALL DLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I), 360: $ X12(I+1,I), LDX12, WORK ) 361: IF( M-P-Q .GE. 1 ) 362: $ CALL DLARF( 'R', M-P-Q, M-Q-I+1, X12(I,I), LDX12, 363: $ TAUQ2(I), X22(Q+1,I), LDX22, WORK ) 364: * 365: END DO 366: * 367: * Reduce columns P + 1, ..., M - Q of X12, X22 368: * 369: DO I = 1, M - P - Q 370: * 371: CALL DSCAL( M-P-Q-I+1, Z2*Z4, X22(Q+I,P+I), LDX22 ) 372: CALL DLARFGP( M-P-Q-I+1, X22(Q+I,P+I), X22(Q+I,P+I+1), 373: $ LDX22, TAUQ2(P+I) ) 374: X22(Q+I,P+I) = ONE 375: CALL DLARF( 'R', M-P-Q-I, M-P-Q-I+1, X22(Q+I,P+I), LDX22, 376: $ TAUQ2(P+I), X22(Q+I+1,P+I), LDX22, WORK ) 377: * 378: END DO 379: * 380: ELSE 381: * 382: * Reduce columns 1, ..., Q of X11, X12, X21, X22 383: * 384: DO I = 1, Q 385: * 386: IF( I .EQ. 1 ) THEN 387: CALL DSCAL( P-I+1, Z1, X11(I,I), LDX11 ) 388: ELSE 389: CALL DSCAL( P-I+1, Z1*COS(PHI(I-1)), X11(I,I), LDX11 ) 390: CALL DAXPY( P-I+1, -Z1*Z3*Z4*SIN(PHI(I-1)), X12(I-1,I), 391: $ LDX12, X11(I,I), LDX11 ) 392: END IF 393: IF( I .EQ. 1 ) THEN 394: CALL DSCAL( M-P-I+1, Z2, X21(I,I), LDX21 ) 395: ELSE 396: CALL DSCAL( M-P-I+1, Z2*COS(PHI(I-1)), X21(I,I), LDX21 ) 397: CALL DAXPY( M-P-I+1, -Z2*Z3*Z4*SIN(PHI(I-1)), X22(I-1,I), 398: $ LDX22, X21(I,I), LDX21 ) 399: END IF 400: * 401: THETA(I) = ATAN2( DNRM2( M-P-I+1, X21(I,I), LDX21 ), 402: $ DNRM2( P-I+1, X11(I,I), LDX11 ) ) 403: * 404: CALL DLARFGP( P-I+1, X11(I,I), X11(I,I+1), LDX11, TAUP1(I) ) 405: X11(I,I) = ONE 406: CALL DLARFGP( M-P-I+1, X21(I,I), X21(I,I+1), LDX21, 407: $ TAUP2(I) ) 408: X21(I,I) = ONE 409: * 410: CALL DLARF( 'R', Q-I, P-I+1, X11(I,I), LDX11, TAUP1(I), 411: $ X11(I+1,I), LDX11, WORK ) 412: CALL DLARF( 'R', M-Q-I+1, P-I+1, X11(I,I), LDX11, TAUP1(I), 413: $ X12(I,I), LDX12, WORK ) 414: CALL DLARF( 'R', Q-I, M-P-I+1, X21(I,I), LDX21, TAUP2(I), 415: $ X21(I+1,I), LDX21, WORK ) 416: CALL DLARF( 'R', M-Q-I+1, M-P-I+1, X21(I,I), LDX21, 417: $ TAUP2(I), X22(I,I), LDX22, WORK ) 418: * 419: IF( I .LT. Q ) THEN 420: CALL DSCAL( Q-I, -Z1*Z3*SIN(THETA(I)), X11(I+1,I), 1 ) 421: CALL DAXPY( Q-I, Z2*Z3*COS(THETA(I)), X21(I+1,I), 1, 422: $ X11(I+1,I), 1 ) 423: END IF 424: CALL DSCAL( M-Q-I+1, -Z1*Z4*SIN(THETA(I)), X12(I,I), 1 ) 425: CALL DAXPY( M-Q-I+1, Z2*Z4*COS(THETA(I)), X22(I,I), 1, 426: $ X12(I,I), 1 ) 427: * 428: IF( I .LT. Q ) 429: $ PHI(I) = ATAN2( DNRM2( Q-I, X11(I+1,I), 1 ), 430: $ DNRM2( M-Q-I+1, X12(I,I), 1 ) ) 431: * 432: IF( I .LT. Q ) THEN 433: CALL DLARFGP( Q-I, X11(I+1,I), X11(I+2,I), 1, TAUQ1(I) ) 434: X11(I+1,I) = ONE 435: END IF 436: CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I+1,I), 1, TAUQ2(I) ) 437: X12(I,I) = ONE 438: * 439: IF( I .LT. Q ) THEN 440: CALL DLARF( 'L', Q-I, P-I, X11(I+1,I), 1, TAUQ1(I), 441: $ X11(I+1,I+1), LDX11, WORK ) 442: CALL DLARF( 'L', Q-I, M-P-I, X11(I+1,I), 1, TAUQ1(I), 443: $ X21(I+1,I+1), LDX21, WORK ) 444: END IF 445: CALL DLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1, TAUQ2(I), 446: $ X12(I,I+1), LDX12, WORK ) 447: CALL DLARF( 'L', M-Q-I+1, M-P-I, X12(I,I), 1, TAUQ2(I), 448: $ X22(I,I+1), LDX22, WORK ) 449: * 450: END DO 451: * 452: * Reduce columns Q + 1, ..., P of X12, X22 453: * 454: DO I = Q + 1, P 455: * 456: CALL DSCAL( M-Q-I+1, -Z1*Z4, X12(I,I), 1 ) 457: CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I+1,I), 1, TAUQ2(I) ) 458: X12(I,I) = ONE 459: * 460: CALL DLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1, TAUQ2(I), 461: $ X12(I,I+1), LDX12, WORK ) 462: IF( M-P-Q .GE. 1 ) 463: $ CALL DLARF( 'L', M-Q-I+1, M-P-Q, X12(I,I), 1, TAUQ2(I), 464: $ X22(I,Q+1), LDX22, WORK ) 465: * 466: END DO 467: * 468: * Reduce columns P + 1, ..., M - Q of X12, X22 469: * 470: DO I = 1, M - P - Q 471: * 472: CALL DSCAL( M-P-Q-I+1, Z2*Z4, X22(P+I,Q+I), 1 ) 473: CALL DLARFGP( M-P-Q-I+1, X22(P+I,Q+I), X22(P+I+1,Q+I), 1, 474: $ TAUQ2(P+I) ) 475: X22(P+I,Q+I) = ONE 476: * 477: CALL DLARF( 'L', M-P-Q-I+1, M-P-Q-I, X22(P+I,Q+I), 1, 478: $ TAUQ2(P+I), X22(P+I,Q+I+1), LDX22, WORK ) 479: * 480: END DO 481: * 482: END IF 483: * 484: RETURN 485: * 486: * End of DORBDB 487: * 488: END 489: