Annotation of rpl/lapack/lapack/zlabrd.f, revision 1.8
1.1 bertrand 1: SUBROUTINE ZLABRD( M, N, NB, A, LDA, D, E, TAUQ, TAUP, X, LDX, Y,
2: $ LDY )
3: *
1.8 ! bertrand 4: * -- LAPACK auxiliary routine (version 3.3.1) --
1.1 bertrand 5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.8 ! bertrand 7: * -- April 2011 --
1.1 bertrand 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
1.8 ! bertrand 23: * transformation Q**H * A * P, and returns the matrices X and Y which
1.1 bertrand 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: *
1.8 ! bertrand 104: * H(i) = I - tauq * v * v**H and G(i) = I - taup * u * u**H
1.1 bertrand 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
1.8 ! bertrand 118: * V and the nb-by-n matrix U**H which are needed, with X and Y, to apply
1.1 bertrand 119: * the transformation to the unreduced part of the matrix, using a block
1.8 ! bertrand 120: * update of the form: A := A - V*Y**H - X*U**H.
1.1 bertrand 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