Annotation of rpl/lapack/lapack/zlabrd.f, revision 1.1
1.1 ! bertrand 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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlabrd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack