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