Annotation of rpl/lapack/lapack/dlabrd.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DLABRD( 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 A( LDA, * ), D( * ), E( * ), TAUP( * ),
! 14: $ TAUQ( * ), X( LDX, * ), Y( LDY, * )
! 15: * ..
! 16: *
! 17: * Purpose
! 18: * =======
! 19: *
! 20: * DLABRD reduces the first NB rows and columns of a real general
! 21: * m by n matrix A to upper or lower bidiagonal form by an orthogonal
! 22: * transformation Q' * A * P, and returns the matrices X and Y which
! 23: * are needed to apply the transformation to the unreduced part of A.
! 24: *
! 25: * If m >= n, A is reduced to upper bidiagonal form; if m < n, to lower
! 26: * bidiagonal form.
! 27: *
! 28: * This is an auxiliary routine called by DGEBRD
! 29: *
! 30: * Arguments
! 31: * =========
! 32: *
! 33: * M (input) INTEGER
! 34: * The number of rows in the matrix A.
! 35: *
! 36: * N (input) INTEGER
! 37: * The number of columns in the matrix A.
! 38: *
! 39: * NB (input) INTEGER
! 40: * The number of leading rows and columns of A to be reduced.
! 41: *
! 42: * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
! 43: * On entry, the m by n general matrix to be reduced.
! 44: * On exit, the first NB rows and columns of the matrix are
! 45: * overwritten; the rest of the array is unchanged.
! 46: * If m >= n, elements on and below the diagonal in the first NB
! 47: * columns, with the array TAUQ, represent the orthogonal
! 48: * matrix Q as a product of elementary reflectors; and
! 49: * elements above the diagonal in the first NB rows, with the
! 50: * array TAUP, represent the orthogonal matrix P as a product
! 51: * of elementary reflectors.
! 52: * If m < n, elements below the diagonal in the first NB
! 53: * columns, with the array TAUQ, represent the orthogonal
! 54: * matrix Q as a product of elementary reflectors, and
! 55: * elements on and above the diagonal in the first NB rows,
! 56: * with the array TAUP, represent the orthogonal matrix P as
! 57: * a product of elementary reflectors.
! 58: * See Further Details.
! 59: *
! 60: * LDA (input) INTEGER
! 61: * The leading dimension of the array A. LDA >= max(1,M).
! 62: *
! 63: * D (output) DOUBLE PRECISION array, dimension (NB)
! 64: * The diagonal elements of the first NB rows and columns of
! 65: * the reduced matrix. D(i) = A(i,i).
! 66: *
! 67: * E (output) DOUBLE PRECISION array, dimension (NB)
! 68: * The off-diagonal elements of the first NB rows and columns of
! 69: * the reduced matrix.
! 70: *
! 71: * TAUQ (output) DOUBLE PRECISION array dimension (NB)
! 72: * The scalar factors of the elementary reflectors which
! 73: * represent the orthogonal matrix Q. See Further Details.
! 74: *
! 75: * TAUP (output) DOUBLE PRECISION array, dimension (NB)
! 76: * The scalar factors of the elementary reflectors which
! 77: * represent the orthogonal matrix P. See Further Details.
! 78: *
! 79: * X (output) DOUBLE PRECISION array, dimension (LDX,NB)
! 80: * The m-by-nb matrix X required to update the unreduced part
! 81: * of A.
! 82: *
! 83: * LDX (input) INTEGER
! 84: * The leading dimension of the array X. LDX >= M.
! 85: *
! 86: * Y (output) DOUBLE PRECISION array, dimension (LDY,NB)
! 87: * The n-by-nb matrix Y required to update the unreduced part
! 88: * of A.
! 89: *
! 90: * LDY (input) INTEGER
! 91: * The leading dimension of the array Y. LDY >= N.
! 92: *
! 93: * Further Details
! 94: * ===============
! 95: *
! 96: * The matrices Q and P are represented as products of elementary
! 97: * reflectors:
! 98: *
! 99: * Q = H(1) H(2) . . . H(nb) and P = G(1) G(2) . . . G(nb)
! 100: *
! 101: * Each H(i) and G(i) has the form:
! 102: *
! 103: * H(i) = I - tauq * v * v' and G(i) = I - taup * u * u'
! 104: *
! 105: * where tauq and taup are real scalars, and v and u are real vectors.
! 106: *
! 107: * If m >= n, v(1:i-1) = 0, v(i) = 1, and v(i:m) is stored on exit in
! 108: * A(i:m,i); u(1:i) = 0, u(i+1) = 1, and u(i+1:n) is stored on exit in
! 109: * A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i).
! 110: *
! 111: * If m < n, v(1:i) = 0, v(i+1) = 1, and v(i+1:m) is stored on exit in
! 112: * A(i+2:m,i); u(1:i-1) = 0, u(i) = 1, and u(i:n) is stored on exit in
! 113: * A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i).
! 114: *
! 115: * The elements of the vectors v and u together form the m-by-nb matrix
! 116: * V and the nb-by-n matrix U' which are needed, with X and Y, to apply
! 117: * the transformation to the unreduced part of the matrix, using a block
! 118: * update of the form: A := A - V*Y' - X*U'.
! 119: *
! 120: * The contents of A on exit are illustrated by the following examples
! 121: * with nb = 2:
! 122: *
! 123: * m = 6 and n = 5 (m > n): m = 5 and n = 6 (m < n):
! 124: *
! 125: * ( 1 1 u1 u1 u1 ) ( 1 u1 u1 u1 u1 u1 )
! 126: * ( v1 1 1 u2 u2 ) ( 1 1 u2 u2 u2 u2 )
! 127: * ( v1 v2 a a a ) ( v1 1 a a a a )
! 128: * ( v1 v2 a a a ) ( v1 v2 a a a a )
! 129: * ( v1 v2 a a a ) ( v1 v2 a a a a )
! 130: * ( v1 v2 a a a )
! 131: *
! 132: * where a denotes an element of the original matrix which is unchanged,
! 133: * vi denotes an element of the vector defining H(i), and ui an element
! 134: * of the vector defining G(i).
! 135: *
! 136: * =====================================================================
! 137: *
! 138: * .. Parameters ..
! 139: DOUBLE PRECISION ZERO, ONE
! 140: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
! 141: * ..
! 142: * .. Local Scalars ..
! 143: INTEGER I
! 144: * ..
! 145: * .. External Subroutines ..
! 146: EXTERNAL DGEMV, DLARFG, DSCAL
! 147: * ..
! 148: * .. Intrinsic Functions ..
! 149: INTRINSIC MIN
! 150: * ..
! 151: * .. Executable Statements ..
! 152: *
! 153: * Quick return if possible
! 154: *
! 155: IF( M.LE.0 .OR. N.LE.0 )
! 156: $ RETURN
! 157: *
! 158: IF( M.GE.N ) THEN
! 159: *
! 160: * Reduce to upper bidiagonal form
! 161: *
! 162: DO 10 I = 1, NB
! 163: *
! 164: * Update A(i:m,i)
! 165: *
! 166: CALL DGEMV( 'No transpose', M-I+1, I-1, -ONE, A( I, 1 ),
! 167: $ LDA, Y( I, 1 ), LDY, ONE, A( I, I ), 1 )
! 168: CALL DGEMV( 'No transpose', M-I+1, I-1, -ONE, X( I, 1 ),
! 169: $ LDX, A( 1, I ), 1, ONE, A( I, I ), 1 )
! 170: *
! 171: * Generate reflection Q(i) to annihilate A(i+1:m,i)
! 172: *
! 173: CALL DLARFG( M-I+1, A( I, I ), A( MIN( I+1, M ), I ), 1,
! 174: $ TAUQ( I ) )
! 175: D( I ) = A( I, I )
! 176: IF( I.LT.N ) THEN
! 177: A( I, I ) = ONE
! 178: *
! 179: * Compute Y(i+1:n,i)
! 180: *
! 181: CALL DGEMV( 'Transpose', M-I+1, N-I, ONE, A( I, I+1 ),
! 182: $ LDA, A( I, I ), 1, ZERO, Y( I+1, I ), 1 )
! 183: CALL DGEMV( 'Transpose', M-I+1, I-1, ONE, A( I, 1 ), LDA,
! 184: $ A( I, I ), 1, ZERO, Y( 1, I ), 1 )
! 185: CALL DGEMV( 'No transpose', N-I, I-1, -ONE, Y( I+1, 1 ),
! 186: $ LDY, Y( 1, I ), 1, ONE, Y( I+1, I ), 1 )
! 187: CALL DGEMV( 'Transpose', M-I+1, I-1, ONE, X( I, 1 ), LDX,
! 188: $ A( I, I ), 1, ZERO, Y( 1, I ), 1 )
! 189: CALL DGEMV( 'Transpose', I-1, N-I, -ONE, A( 1, I+1 ),
! 190: $ LDA, Y( 1, I ), 1, ONE, Y( I+1, I ), 1 )
! 191: CALL DSCAL( N-I, TAUQ( I ), Y( I+1, I ), 1 )
! 192: *
! 193: * Update A(i,i+1:n)
! 194: *
! 195: CALL DGEMV( 'No transpose', N-I, I, -ONE, Y( I+1, 1 ),
! 196: $ LDY, A( I, 1 ), LDA, ONE, A( I, I+1 ), LDA )
! 197: CALL DGEMV( 'Transpose', I-1, N-I, -ONE, A( 1, I+1 ),
! 198: $ LDA, X( I, 1 ), LDX, ONE, A( I, I+1 ), LDA )
! 199: *
! 200: * Generate reflection P(i) to annihilate A(i,i+2:n)
! 201: *
! 202: CALL DLARFG( N-I, A( I, I+1 ), A( I, MIN( I+2, N ) ),
! 203: $ LDA, TAUP( I ) )
! 204: E( I ) = A( I, I+1 )
! 205: A( I, I+1 ) = ONE
! 206: *
! 207: * Compute X(i+1:m,i)
! 208: *
! 209: CALL DGEMV( 'No transpose', M-I, N-I, ONE, A( I+1, I+1 ),
! 210: $ LDA, A( I, I+1 ), LDA, ZERO, X( I+1, I ), 1 )
! 211: CALL DGEMV( 'Transpose', N-I, I, ONE, Y( I+1, 1 ), LDY,
! 212: $ A( I, I+1 ), LDA, ZERO, X( 1, I ), 1 )
! 213: CALL DGEMV( 'No transpose', M-I, I, -ONE, A( I+1, 1 ),
! 214: $ LDA, X( 1, I ), 1, ONE, X( I+1, I ), 1 )
! 215: CALL DGEMV( 'No transpose', I-1, N-I, ONE, A( 1, I+1 ),
! 216: $ LDA, A( I, I+1 ), LDA, ZERO, X( 1, I ), 1 )
! 217: CALL DGEMV( 'No transpose', M-I, I-1, -ONE, X( I+1, 1 ),
! 218: $ LDX, X( 1, I ), 1, ONE, X( I+1, I ), 1 )
! 219: CALL DSCAL( M-I, TAUP( I ), X( I+1, I ), 1 )
! 220: END IF
! 221: 10 CONTINUE
! 222: ELSE
! 223: *
! 224: * Reduce to lower bidiagonal form
! 225: *
! 226: DO 20 I = 1, NB
! 227: *
! 228: * Update A(i,i:n)
! 229: *
! 230: CALL DGEMV( 'No transpose', N-I+1, I-1, -ONE, Y( I, 1 ),
! 231: $ LDY, A( I, 1 ), LDA, ONE, A( I, I ), LDA )
! 232: CALL DGEMV( 'Transpose', I-1, N-I+1, -ONE, A( 1, I ), LDA,
! 233: $ X( I, 1 ), LDX, ONE, A( I, I ), LDA )
! 234: *
! 235: * Generate reflection P(i) to annihilate A(i,i+1:n)
! 236: *
! 237: CALL DLARFG( N-I+1, A( I, I ), A( I, MIN( I+1, N ) ), LDA,
! 238: $ TAUP( I ) )
! 239: D( I ) = A( I, I )
! 240: IF( I.LT.M ) THEN
! 241: A( I, I ) = ONE
! 242: *
! 243: * Compute X(i+1:m,i)
! 244: *
! 245: CALL DGEMV( 'No transpose', M-I, N-I+1, ONE, A( I+1, I ),
! 246: $ LDA, A( I, I ), LDA, ZERO, X( I+1, I ), 1 )
! 247: CALL DGEMV( 'Transpose', N-I+1, I-1, ONE, Y( I, 1 ), LDY,
! 248: $ A( I, I ), LDA, ZERO, X( 1, I ), 1 )
! 249: CALL DGEMV( 'No transpose', M-I, I-1, -ONE, A( I+1, 1 ),
! 250: $ LDA, X( 1, I ), 1, ONE, X( I+1, I ), 1 )
! 251: CALL DGEMV( 'No transpose', I-1, N-I+1, ONE, A( 1, I ),
! 252: $ LDA, A( I, I ), LDA, ZERO, X( 1, I ), 1 )
! 253: CALL DGEMV( 'No transpose', M-I, I-1, -ONE, X( I+1, 1 ),
! 254: $ LDX, X( 1, I ), 1, ONE, X( I+1, I ), 1 )
! 255: CALL DSCAL( M-I, TAUP( I ), X( I+1, I ), 1 )
! 256: *
! 257: * Update A(i+1:m,i)
! 258: *
! 259: CALL DGEMV( 'No transpose', M-I, I-1, -ONE, A( I+1, 1 ),
! 260: $ LDA, Y( I, 1 ), LDY, ONE, A( I+1, I ), 1 )
! 261: CALL DGEMV( 'No transpose', M-I, I, -ONE, X( I+1, 1 ),
! 262: $ LDX, A( 1, I ), 1, ONE, A( I+1, I ), 1 )
! 263: *
! 264: * Generate reflection Q(i) to annihilate A(i+2:m,i)
! 265: *
! 266: CALL DLARFG( M-I, A( I+1, I ), A( MIN( I+2, M ), I ), 1,
! 267: $ TAUQ( I ) )
! 268: E( I ) = A( I+1, I )
! 269: A( I+1, I ) = ONE
! 270: *
! 271: * Compute Y(i+1:n,i)
! 272: *
! 273: CALL DGEMV( 'Transpose', M-I, N-I, ONE, A( I+1, I+1 ),
! 274: $ LDA, A( I+1, I ), 1, ZERO, Y( I+1, I ), 1 )
! 275: CALL DGEMV( 'Transpose', M-I, I-1, ONE, A( I+1, 1 ), LDA,
! 276: $ A( I+1, I ), 1, ZERO, Y( 1, I ), 1 )
! 277: CALL DGEMV( 'No transpose', N-I, I-1, -ONE, Y( I+1, 1 ),
! 278: $ LDY, Y( 1, I ), 1, ONE, Y( I+1, I ), 1 )
! 279: CALL DGEMV( 'Transpose', M-I, I, ONE, X( I+1, 1 ), LDX,
! 280: $ A( I+1, I ), 1, ZERO, Y( 1, I ), 1 )
! 281: CALL DGEMV( 'Transpose', I, N-I, -ONE, A( 1, I+1 ), LDA,
! 282: $ Y( 1, I ), 1, ONE, Y( I+1, I ), 1 )
! 283: CALL DSCAL( N-I, TAUQ( I ), Y( I+1, I ), 1 )
! 284: END IF
! 285: 20 CONTINUE
! 286: END IF
! 287: RETURN
! 288: *
! 289: * End of DLABRD
! 290: *
! 291: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlabrd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack