001:       SUBROUTINE DLABRD( M, N, NB, A, LDA, D, E, TAUQ, TAUP, X, LDX, Y,
002:      $                   LDY )
003: *
004: *  -- LAPACK auxiliary routine (version 3.2) --
005: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       INTEGER            LDA, LDX, LDY, M, N, NB
011: *     ..
012: *     .. Array Arguments ..
013:       DOUBLE PRECISION   A( LDA, * ), D( * ), E( * ), TAUP( * ),
014:      $                   TAUQ( * ), X( LDX, * ), Y( LDY, * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  DLABRD reduces the first NB rows and columns of a real general
021: *  m by n matrix A to upper or lower bidiagonal form by an orthogonal
022: *  transformation Q' * A * P, and returns the matrices X and Y which
023: *  are needed to apply the transformation to the unreduced part of A.
024: *
025: *  If m >= n, A is reduced to upper bidiagonal form; if m < n, to lower
026: *  bidiagonal form.
027: *
028: *  This is an auxiliary routine called by DGEBRD
029: *
030: *  Arguments
031: *  =========
032: *
033: *  M       (input) INTEGER
034: *          The number of rows in the matrix A.
035: *
036: *  N       (input) INTEGER
037: *          The number of columns in the matrix A.
038: *
039: *  NB      (input) INTEGER
040: *          The number of leading rows and columns of A to be reduced.
041: *
042: *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
043: *          On entry, the m by n general matrix to be reduced.
044: *          On exit, the first NB rows and columns of the matrix are
045: *          overwritten; the rest of the array is unchanged.
046: *          If m >= n, elements on and below the diagonal in the first NB
047: *            columns, with the array TAUQ, represent the orthogonal
048: *            matrix Q as a product of elementary reflectors; and
049: *            elements above the diagonal in the first NB rows, with the
050: *            array TAUP, represent the orthogonal matrix P as a product
051: *            of elementary reflectors.
052: *          If m < n, elements below the diagonal in the first NB
053: *            columns, with the array TAUQ, represent the orthogonal
054: *            matrix Q as a product of elementary reflectors, and
055: *            elements on and above the diagonal in the first NB rows,
056: *            with the array TAUP, represent the orthogonal matrix P as
057: *            a product of elementary reflectors.
058: *          See Further Details.
059: *
060: *  LDA     (input) INTEGER
061: *          The leading dimension of the array A.  LDA >= max(1,M).
062: *
063: *  D       (output) DOUBLE PRECISION array, dimension (NB)
064: *          The diagonal elements of the first NB rows and columns of
065: *          the reduced matrix.  D(i) = A(i,i).
066: *
067: *  E       (output) DOUBLE PRECISION array, dimension (NB)
068: *          The off-diagonal elements of the first NB rows and columns of
069: *          the reduced matrix.
070: *
071: *  TAUQ    (output) DOUBLE PRECISION array dimension (NB)
072: *          The scalar factors of the elementary reflectors which
073: *          represent the orthogonal matrix Q. See Further Details.
074: *
075: *  TAUP    (output) DOUBLE PRECISION array, dimension (NB)
076: *          The scalar factors of the elementary reflectors which
077: *          represent the orthogonal matrix P. See Further Details.
078: *
079: *  X       (output) DOUBLE PRECISION array, dimension (LDX,NB)
080: *          The m-by-nb matrix X required to update the unreduced part
081: *          of A.
082: *
083: *  LDX     (input) INTEGER
084: *          The leading dimension of the array X. LDX >= M.
085: *
086: *  Y       (output) DOUBLE PRECISION array, dimension (LDY,NB)
087: *          The n-by-nb matrix Y required to update the unreduced part
088: *          of A.
089: *
090: *  LDY     (input) INTEGER
091: *          The leading dimension of the array Y. LDY >= N.
092: *
093: *  Further Details
094: *  ===============
095: *
096: *  The matrices Q and P are represented as products of elementary
097: *  reflectors:
098: *
099: *     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
292: