![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZLABRD( M, N, NB, A, LDA, D, E, TAUQ, TAUP, X, LDX, Y, 2: $ LDY ) 3: * 4: * -- LAPACK auxiliary routine (version 3.2) -- 5: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 7: * November 2006 8: * 9: * .. Scalar Arguments .. 10: INTEGER LDA, LDX, LDY, M, N, NB 11: * .. 12: * .. Array Arguments .. 13: DOUBLE PRECISION D( * ), E( * ) 14: COMPLEX*16 A( LDA, * ), TAUP( * ), TAUQ( * ), X( LDX, * ), 15: $ Y( LDY, * ) 16: * .. 17: * 18: * Purpose 19: * ======= 20: * 21: * ZLABRD reduces the first NB rows and columns of a complex general 22: * m by n matrix A to upper or lower real bidiagonal form by a unitary 23: * transformation Q' * A * P, and returns the matrices X and Y which 24: * are needed to apply the transformation to the unreduced part of A. 25: * 26: * If m >= n, A is reduced to upper bidiagonal form; if m < n, to lower 27: * bidiagonal form. 28: * 29: * This is an auxiliary routine called by ZGEBRD 30: * 31: * Arguments 32: * ========= 33: * 34: * M (input) INTEGER 35: * The number of rows in the matrix A. 36: * 37: * N (input) INTEGER 38: * The number of columns in the matrix A. 39: * 40: * NB (input) INTEGER 41: * The number of leading rows and columns of A to be reduced. 42: * 43: * A (input/output) COMPLEX*16 array, dimension (LDA,N) 44: * On entry, the m by n general matrix to be reduced. 45: * On exit, the first NB rows and columns of the matrix are 46: * overwritten; the rest of the array is unchanged. 47: * If m >= n, elements on and below the diagonal in the first NB 48: * columns, with the array TAUQ, represent the unitary 49: * matrix Q as a product of elementary reflectors; and 50: * elements above the diagonal in the first NB rows, with the 51: * array TAUP, represent the unitary matrix P as a product 52: * of elementary reflectors. 53: * If m < n, elements below the diagonal in the first NB 54: * columns, with the array TAUQ, represent the unitary 55: * matrix Q as a product of elementary reflectors, and 56: * elements on and above the diagonal in the first NB rows, 57: * with the array TAUP, represent the unitary matrix P as 58: * a product of elementary reflectors. 59: * See Further Details. 60: * 61: * LDA (input) INTEGER 62: * The leading dimension of the array A. LDA >= max(1,M). 63: * 64: * D (output) DOUBLE PRECISION array, dimension (NB) 65: * The diagonal elements of the first NB rows and columns of 66: * the reduced matrix. D(i) = A(i,i). 67: * 68: * E (output) DOUBLE PRECISION array, dimension (NB) 69: * The off-diagonal elements of the first NB rows and columns of 70: * the reduced matrix. 71: * 72: * TAUQ (output) COMPLEX*16 array dimension (NB) 73: * The scalar factors of the elementary reflectors which 74: * represent the unitary matrix Q. See Further Details. 75: * 76: * TAUP (output) COMPLEX*16 array, dimension (NB) 77: * The scalar factors of the elementary reflectors which 78: * represent the unitary matrix P. See Further Details. 79: * 80: * X (output) COMPLEX*16 array, dimension (LDX,NB) 81: * The m-by-nb matrix X required to update the unreduced part 82: * of A. 83: * 84: * LDX (input) INTEGER 85: * The leading dimension of the array X. LDX >= max(1,M). 86: * 87: * Y (output) COMPLEX*16 array, dimension (LDY,NB) 88: * The n-by-nb matrix Y required to update the unreduced part 89: * of A. 90: * 91: * LDY (input) INTEGER 92: * The leading dimension of the array Y. LDY >= max(1,N). 93: * 94: * Further Details 95: * =============== 96: * 97: * The matrices Q and P are represented as products of elementary 98: * reflectors: 99: * 100: * Q = H(1) H(2) . . . H(nb) and P = G(1) G(2) . . . G(nb) 101: * 102: * Each H(i) and G(i) has the form: 103: * 104: * H(i) = I - tauq * v * v' and G(i) = I - taup * u * u' 105: * 106: * where tauq and taup are complex scalars, and v and u are complex 107: * vectors. 108: * 109: * If m >= n, v(1:i-1) = 0, v(i) = 1, and v(i:m) is stored on exit in 110: * A(i:m,i); u(1:i) = 0, u(i+1) = 1, and u(i+1:n) is stored on exit in 111: * A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i). 112: * 113: * If m < n, v(1:i) = 0, v(i+1) = 1, and v(i+1:m) is stored on exit in 114: * A(i+2:m,i); u(1:i-1) = 0, u(i) = 1, and u(i:n) is stored on exit in 115: * A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i). 116: * 117: * The elements of the vectors v and u together form the m-by-nb matrix 118: * V and the nb-by-n matrix U' which are needed, with X and Y, to apply 119: * the transformation to the unreduced part of the matrix, using a block 120: * update of the form: A := A - V*Y' - X*U'. 121: * 122: * The contents of A on exit are illustrated by the following examples 123: * with nb = 2: 124: * 125: * m = 6 and n = 5 (m > n): m = 5 and n = 6 (m < n): 126: * 127: * ( 1 1 u1 u1 u1 ) ( 1 u1 u1 u1 u1 u1 ) 128: * ( v1 1 1 u2 u2 ) ( 1 1 u2 u2 u2 u2 ) 129: * ( v1 v2 a a a ) ( v1 1 a a a a ) 130: * ( v1 v2 a a a ) ( v1 v2 a a a a ) 131: * ( v1 v2 a a a ) ( v1 v2 a a a a ) 132: * ( v1 v2 a a a ) 133: * 134: * where a denotes an element of the original matrix which is unchanged, 135: * vi denotes an element of the vector defining H(i), and ui an element 136: * of the vector defining G(i). 137: * 138: * ===================================================================== 139: * 140: * .. Parameters .. 141: COMPLEX*16 ZERO, ONE 142: PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ), 143: $ ONE = ( 1.0D+0, 0.0D+0 ) ) 144: * .. 145: * .. Local Scalars .. 146: INTEGER I 147: COMPLEX*16 ALPHA 148: * .. 149: * .. External Subroutines .. 150: EXTERNAL ZGEMV, ZLACGV, ZLARFG, ZSCAL 151: * .. 152: * .. Intrinsic Functions .. 153: INTRINSIC MIN 154: * .. 155: * .. Executable Statements .. 156: * 157: * Quick return if possible 158: * 159: IF( M.LE.0 .OR. N.LE.0 ) 160: $ RETURN 161: * 162: IF( M.GE.N ) THEN 163: * 164: * Reduce to upper bidiagonal form 165: * 166: DO 10 I = 1, NB 167: * 168: * Update A(i:m,i) 169: * 170: CALL ZLACGV( I-1, Y( I, 1 ), LDY ) 171: CALL ZGEMV( 'No transpose', M-I+1, I-1, -ONE, A( I, 1 ), 172: $ LDA, Y( I, 1 ), LDY, ONE, A( I, I ), 1 ) 173: CALL ZLACGV( I-1, Y( I, 1 ), LDY ) 174: CALL ZGEMV( 'No transpose', M-I+1, I-1, -ONE, X( I, 1 ), 175: $ LDX, A( 1, I ), 1, ONE, A( I, I ), 1 ) 176: * 177: * Generate reflection Q(i) to annihilate A(i+1:m,i) 178: * 179: ALPHA = A( I, I ) 180: CALL ZLARFG( M-I+1, ALPHA, A( MIN( I+1, M ), I ), 1, 181: $ TAUQ( I ) ) 182: D( I ) = ALPHA 183: IF( I.LT.N ) THEN 184: A( I, I ) = ONE 185: * 186: * Compute Y(i+1:n,i) 187: * 188: CALL ZGEMV( 'Conjugate transpose', M-I+1, N-I, ONE, 189: $ A( I, I+1 ), LDA, A( I, I ), 1, ZERO, 190: $ Y( I+1, I ), 1 ) 191: CALL ZGEMV( 'Conjugate transpose', M-I+1, I-1, ONE, 192: $ A( I, 1 ), LDA, A( I, I ), 1, ZERO, 193: $ Y( 1, I ), 1 ) 194: CALL ZGEMV( 'No transpose', N-I, I-1, -ONE, Y( I+1, 1 ), 195: $ LDY, Y( 1, I ), 1, ONE, Y( I+1, I ), 1 ) 196: CALL ZGEMV( 'Conjugate transpose', M-I+1, I-1, ONE, 197: $ X( I, 1 ), LDX, A( I, I ), 1, ZERO, 198: $ Y( 1, I ), 1 ) 199: CALL ZGEMV( 'Conjugate transpose', I-1, N-I, -ONE, 200: $ A( 1, I+1 ), LDA, Y( 1, I ), 1, ONE, 201: $ Y( I+1, I ), 1 ) 202: CALL ZSCAL( N-I, TAUQ( I ), Y( I+1, I ), 1 ) 203: * 204: * Update A(i,i+1:n) 205: * 206: CALL ZLACGV( N-I, A( I, I+1 ), LDA ) 207: CALL ZLACGV( I, A( I, 1 ), LDA ) 208: CALL ZGEMV( 'No transpose', N-I, I, -ONE, Y( I+1, 1 ), 209: $ LDY, A( I, 1 ), LDA, ONE, A( I, I+1 ), LDA ) 210: CALL ZLACGV( I, A( I, 1 ), LDA ) 211: CALL ZLACGV( I-1, X( I, 1 ), LDX ) 212: CALL ZGEMV( 'Conjugate transpose', I-1, N-I, -ONE, 213: $ A( 1, I+1 ), LDA, X( I, 1 ), LDX, ONE, 214: $ A( I, I+1 ), LDA ) 215: CALL ZLACGV( I-1, X( I, 1 ), LDX ) 216: * 217: * Generate reflection P(i) to annihilate A(i,i+2:n) 218: * 219: ALPHA = A( I, I+1 ) 220: CALL ZLARFG( N-I, ALPHA, A( I, MIN( I+2, N ) ), LDA, 221: $ TAUP( I ) ) 222: E( I ) = ALPHA 223: A( I, I+1 ) = ONE 224: * 225: * Compute X(i+1:m,i) 226: * 227: CALL ZGEMV( 'No transpose', M-I, N-I, ONE, A( I+1, I+1 ), 228: $ LDA, A( I, I+1 ), LDA, ZERO, X( I+1, I ), 1 ) 229: CALL ZGEMV( 'Conjugate transpose', N-I, I, ONE, 230: $ Y( I+1, 1 ), LDY, A( I, I+1 ), LDA, ZERO, 231: $ X( 1, I ), 1 ) 232: CALL ZGEMV( 'No transpose', M-I, I, -ONE, A( I+1, 1 ), 233: $ LDA, X( 1, I ), 1, ONE, X( I+1, I ), 1 ) 234: CALL ZGEMV( 'No transpose', I-1, N-I, ONE, A( 1, I+1 ), 235: $ LDA, A( I, I+1 ), LDA, ZERO, X( 1, I ), 1 ) 236: CALL ZGEMV( 'No transpose', M-I, I-1, -ONE, X( I+1, 1 ), 237: $ LDX, X( 1, I ), 1, ONE, X( I+1, I ), 1 ) 238: CALL ZSCAL( M-I, TAUP( I ), X( I+1, I ), 1 ) 239: CALL ZLACGV( N-I, A( I, I+1 ), LDA ) 240: END IF 241: 10 CONTINUE 242: ELSE 243: * 244: * Reduce to lower bidiagonal form 245: * 246: DO 20 I = 1, NB 247: * 248: * Update A(i,i:n) 249: * 250: CALL ZLACGV( N-I+1, A( I, I ), LDA ) 251: CALL ZLACGV( I-1, A( I, 1 ), LDA ) 252: CALL ZGEMV( 'No transpose', N-I+1, I-1, -ONE, Y( I, 1 ), 253: $ LDY, A( I, 1 ), LDA, ONE, A( I, I ), LDA ) 254: CALL ZLACGV( I-1, A( I, 1 ), LDA ) 255: CALL ZLACGV( I-1, X( I, 1 ), LDX ) 256: CALL ZGEMV( 'Conjugate transpose', I-1, N-I+1, -ONE, 257: $ A( 1, I ), LDA, X( I, 1 ), LDX, ONE, A( I, I ), 258: $ LDA ) 259: CALL ZLACGV( I-1, X( I, 1 ), LDX ) 260: * 261: * Generate reflection P(i) to annihilate A(i,i+1:n) 262: * 263: ALPHA = A( I, I ) 264: CALL ZLARFG( N-I+1, ALPHA, A( I, MIN( I+1, N ) ), LDA, 265: $ TAUP( I ) ) 266: D( I ) = ALPHA 267: IF( I.LT.M ) THEN 268: A( I, I ) = ONE 269: * 270: * Compute X(i+1:m,i) 271: * 272: CALL ZGEMV( 'No transpose', M-I, N-I+1, ONE, A( I+1, I ), 273: $ LDA, A( I, I ), LDA, ZERO, X( I+1, I ), 1 ) 274: CALL ZGEMV( 'Conjugate transpose', N-I+1, I-1, ONE, 275: $ Y( I, 1 ), LDY, A( I, I ), LDA, ZERO, 276: $ X( 1, I ), 1 ) 277: CALL ZGEMV( 'No transpose', M-I, I-1, -ONE, A( I+1, 1 ), 278: $ LDA, X( 1, I ), 1, ONE, X( I+1, I ), 1 ) 279: CALL ZGEMV( 'No transpose', I-1, N-I+1, ONE, A( 1, I ), 280: $ LDA, A( I, I ), LDA, ZERO, X( 1, I ), 1 ) 281: CALL ZGEMV( 'No transpose', M-I, I-1, -ONE, X( I+1, 1 ), 282: $ LDX, X( 1, I ), 1, ONE, X( I+1, I ), 1 ) 283: CALL ZSCAL( M-I, TAUP( I ), X( I+1, I ), 1 ) 284: CALL ZLACGV( N-I+1, A( I, I ), LDA ) 285: * 286: * Update A(i+1:m,i) 287: * 288: CALL ZLACGV( I-1, Y( I, 1 ), LDY ) 289: CALL ZGEMV( 'No transpose', M-I, I-1, -ONE, A( I+1, 1 ), 290: $ LDA, Y( I, 1 ), LDY, ONE, A( I+1, I ), 1 ) 291: CALL ZLACGV( I-1, Y( I, 1 ), LDY ) 292: CALL ZGEMV( 'No transpose', M-I, I, -ONE, X( I+1, 1 ), 293: $ LDX, A( 1, I ), 1, ONE, A( I+1, I ), 1 ) 294: * 295: * Generate reflection Q(i) to annihilate A(i+2:m,i) 296: * 297: ALPHA = A( I+1, I ) 298: CALL ZLARFG( M-I, ALPHA, A( MIN( I+2, M ), I ), 1, 299: $ TAUQ( I ) ) 300: E( I ) = ALPHA 301: A( I+1, I ) = ONE 302: * 303: * Compute Y(i+1:n,i) 304: * 305: CALL ZGEMV( 'Conjugate transpose', M-I, N-I, ONE, 306: $ A( I+1, I+1 ), LDA, A( I+1, I ), 1, ZERO, 307: $ Y( I+1, I ), 1 ) 308: CALL ZGEMV( 'Conjugate transpose', M-I, I-1, ONE, 309: $ A( I+1, 1 ), LDA, A( I+1, I ), 1, ZERO, 310: $ Y( 1, I ), 1 ) 311: CALL ZGEMV( 'No transpose', N-I, I-1, -ONE, Y( I+1, 1 ), 312: $ LDY, Y( 1, I ), 1, ONE, Y( I+1, I ), 1 ) 313: CALL ZGEMV( 'Conjugate transpose', M-I, I, ONE, 314: $ X( I+1, 1 ), LDX, A( I+1, I ), 1, ZERO, 315: $ Y( 1, I ), 1 ) 316: CALL ZGEMV( 'Conjugate transpose', I, N-I, -ONE, 317: $ A( 1, I+1 ), LDA, Y( 1, I ), 1, ONE, 318: $ Y( I+1, I ), 1 ) 319: CALL ZSCAL( N-I, TAUQ( I ), Y( I+1, I ), 1 ) 320: ELSE 321: CALL ZLACGV( N-I+1, A( I, I ), LDA ) 322: END IF 323: 20 CONTINUE 324: END IF 325: RETURN 326: * 327: * End of ZLABRD 328: * 329: END