1: *> \brief \b ZLAHRD reduces the first nb columns of a general rectangular matrix A so that elements below the k-th subdiagonal are zero, and returns auxiliary matrices which are needed to apply the transformation to the unreduced part of A.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZLAHRD + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlahrd.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlahrd.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlahrd.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZLAHRD( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY )
22: *
23: * .. Scalar Arguments ..
24: * INTEGER K, LDA, LDT, LDY, N, NB
25: * ..
26: * .. Array Arguments ..
27: * COMPLEX*16 A( LDA, * ), T( LDT, NB ), TAU( NB ),
28: * $ Y( LDY, NB )
29: * ..
30: *
31: *
32: *> \par Purpose:
33: * =============
34: *>
35: *> \verbatim
36: *>
37: *> This routine is deprecated and has been replaced by routine ZLAHR2.
38: *>
39: *> ZLAHRD reduces the first NB columns of a complex general n-by-(n-k+1)
40: *> matrix A so that elements below the k-th subdiagonal are zero. The
41: *> reduction is performed by a unitary similarity transformation
42: *> Q**H * A * Q. The routine returns the matrices V and T which determine
43: *> Q as a block reflector I - V*T*V**H, and also the matrix Y = A * V * T.
44: *> \endverbatim
45: *
46: * Arguments:
47: * ==========
48: *
49: *> \param[in] N
50: *> \verbatim
51: *> N is INTEGER
52: *> The order of the matrix A.
53: *> \endverbatim
54: *>
55: *> \param[in] K
56: *> \verbatim
57: *> K is INTEGER
58: *> The offset for the reduction. Elements below the k-th
59: *> subdiagonal in the first NB columns are reduced to zero.
60: *> \endverbatim
61: *>
62: *> \param[in] NB
63: *> \verbatim
64: *> NB is INTEGER
65: *> The number of columns to be reduced.
66: *> \endverbatim
67: *>
68: *> \param[in,out] A
69: *> \verbatim
70: *> A is COMPLEX*16 array, dimension (LDA,N-K+1)
71: *> On entry, the n-by-(n-k+1) general matrix A.
72: *> On exit, the elements on and above the k-th subdiagonal in
73: *> the first NB columns are overwritten with the corresponding
74: *> elements of the reduced matrix; the elements below the k-th
75: *> subdiagonal, with the array TAU, represent the matrix Q as a
76: *> product of elementary reflectors. The other columns of A are
77: *> unchanged. See Further Details.
78: *> \endverbatim
79: *>
80: *> \param[in] LDA
81: *> \verbatim
82: *> LDA is INTEGER
83: *> The leading dimension of the array A. LDA >= max(1,N).
84: *> \endverbatim
85: *>
86: *> \param[out] TAU
87: *> \verbatim
88: *> TAU is COMPLEX*16 array, dimension (NB)
89: *> The scalar factors of the elementary reflectors. See Further
90: *> Details.
91: *> \endverbatim
92: *>
93: *> \param[out] T
94: *> \verbatim
95: *> T is COMPLEX*16 array, dimension (LDT,NB)
96: *> The upper triangular matrix T.
97: *> \endverbatim
98: *>
99: *> \param[in] LDT
100: *> \verbatim
101: *> LDT is INTEGER
102: *> The leading dimension of the array T. LDT >= NB.
103: *> \endverbatim
104: *>
105: *> \param[out] Y
106: *> \verbatim
107: *> Y is COMPLEX*16 array, dimension (LDY,NB)
108: *> The n-by-nb matrix Y.
109: *> \endverbatim
110: *>
111: *> \param[in] LDY
112: *> \verbatim
113: *> LDY is INTEGER
114: *> The leading dimension of the array Y. LDY >= max(1,N).
115: *> \endverbatim
116: *
117: * Authors:
118: * ========
119: *
120: *> \author Univ. of Tennessee
121: *> \author Univ. of California Berkeley
122: *> \author Univ. of Colorado Denver
123: *> \author NAG Ltd.
124: *
125: *> \ingroup complex16OTHERauxiliary
126: *
127: *> \par Further Details:
128: * =====================
129: *>
130: *> \verbatim
131: *>
132: *> The matrix Q is represented as a product of nb elementary reflectors
133: *>
134: *> Q = H(1) H(2) . . . H(nb).
135: *>
136: *> Each H(i) has the form
137: *>
138: *> H(i) = I - tau * v * v**H
139: *>
140: *> where tau is a complex scalar, and v is a complex vector with
141: *> v(1:i+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in
142: *> A(i+k+1:n,i), and tau in TAU(i).
143: *>
144: *> The elements of the vectors v together form the (n-k+1)-by-nb matrix
145: *> V which is needed, with T and Y, to apply the transformation to the
146: *> unreduced part of the matrix, using an update of the form:
147: *> A := (I - V*T*V**H) * (A - Y*V**H).
148: *>
149: *> The contents of A on exit are illustrated by the following example
150: *> with n = 7, k = 3 and nb = 2:
151: *>
152: *> ( a h a a a )
153: *> ( a h a a a )
154: *> ( a h a a a )
155: *> ( h h a a a )
156: *> ( v1 h a a a )
157: *> ( v1 v2 a a a )
158: *> ( v1 v2 a a a )
159: *>
160: *> where a denotes an element of the original matrix A, h denotes a
161: *> modified element of the upper Hessenberg matrix H, and vi denotes an
162: *> element of the vector defining H(i).
163: *> \endverbatim
164: *>
165: * =====================================================================
166: SUBROUTINE ZLAHRD( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY )
167: *
168: * -- LAPACK auxiliary routine --
169: * -- LAPACK is a software package provided by Univ. of Tennessee, --
170: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
171: *
172: * .. Scalar Arguments ..
173: INTEGER K, LDA, LDT, LDY, N, NB
174: * ..
175: * .. Array Arguments ..
176: COMPLEX*16 A( LDA, * ), T( LDT, NB ), TAU( NB ),
177: $ Y( LDY, NB )
178: * ..
179: *
180: * =====================================================================
181: *
182: * .. Parameters ..
183: COMPLEX*16 ZERO, ONE
184: PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ),
185: $ ONE = ( 1.0D+0, 0.0D+0 ) )
186: * ..
187: * .. Local Scalars ..
188: INTEGER I
189: COMPLEX*16 EI
190: * ..
191: * .. External Subroutines ..
192: EXTERNAL ZAXPY, ZCOPY, ZGEMV, ZLACGV, ZLARFG, ZSCAL,
193: $ ZTRMV
194: * ..
195: * .. Intrinsic Functions ..
196: INTRINSIC MIN
197: * ..
198: * .. Executable Statements ..
199: *
200: * Quick return if possible
201: *
202: IF( N.LE.1 )
203: $ RETURN
204: *
205: DO 10 I = 1, NB
206: IF( I.GT.1 ) THEN
207: *
208: * Update A(1:n,i)
209: *
210: * Compute i-th column of A - Y * V**H
211: *
212: CALL ZLACGV( I-1, A( K+I-1, 1 ), LDA )
213: CALL ZGEMV( 'No transpose', N, I-1, -ONE, Y, LDY,
214: $ A( K+I-1, 1 ), LDA, ONE, A( 1, I ), 1 )
215: CALL ZLACGV( I-1, A( K+I-1, 1 ), LDA )
216: *
217: * Apply I - V * T**H * V**H to this column (call it b) from the
218: * left, using the last column of T as workspace
219: *
220: * Let V = ( V1 ) and b = ( b1 ) (first I-1 rows)
221: * ( V2 ) ( b2 )
222: *
223: * where V1 is unit lower triangular
224: *
225: * w := V1**H * b1
226: *
227: CALL ZCOPY( I-1, A( K+1, I ), 1, T( 1, NB ), 1 )
228: CALL ZTRMV( 'Lower', 'Conjugate transpose', 'Unit', I-1,
229: $ A( K+1, 1 ), LDA, T( 1, NB ), 1 )
230: *
231: * w := w + V2**H *b2
232: *
233: CALL ZGEMV( 'Conjugate transpose', N-K-I+1, I-1, ONE,
234: $ A( K+I, 1 ), LDA, A( K+I, I ), 1, ONE,
235: $ T( 1, NB ), 1 )
236: *
237: * w := T**H *w
238: *
239: CALL ZTRMV( 'Upper', 'Conjugate transpose', 'Non-unit', I-1,
240: $ T, LDT, T( 1, NB ), 1 )
241: *
242: * b2 := b2 - V2*w
243: *
244: CALL ZGEMV( 'No transpose', N-K-I+1, I-1, -ONE, A( K+I, 1 ),
245: $ LDA, T( 1, NB ), 1, ONE, A( K+I, I ), 1 )
246: *
247: * b1 := b1 - V1*w
248: *
249: CALL ZTRMV( 'Lower', 'No transpose', 'Unit', I-1,
250: $ A( K+1, 1 ), LDA, T( 1, NB ), 1 )
251: CALL ZAXPY( I-1, -ONE, T( 1, NB ), 1, A( K+1, I ), 1 )
252: *
253: A( K+I-1, I-1 ) = EI
254: END IF
255: *
256: * Generate the elementary reflector H(i) to annihilate
257: * A(k+i+1:n,i)
258: *
259: EI = A( K+I, I )
260: CALL ZLARFG( N-K-I+1, EI, A( MIN( K+I+1, N ), I ), 1,
261: $ TAU( I ) )
262: A( K+I, I ) = ONE
263: *
264: * Compute Y(1:n,i)
265: *
266: CALL ZGEMV( 'No transpose', N, N-K-I+1, ONE, A( 1, I+1 ), LDA,
267: $ A( K+I, I ), 1, ZERO, Y( 1, I ), 1 )
268: CALL ZGEMV( 'Conjugate transpose', N-K-I+1, I-1, ONE,
269: $ A( K+I, 1 ), LDA, A( K+I, I ), 1, ZERO, T( 1, I ),
270: $ 1 )
271: CALL ZGEMV( 'No transpose', N, I-1, -ONE, Y, LDY, T( 1, I ), 1,
272: $ ONE, Y( 1, I ), 1 )
273: CALL ZSCAL( N, TAU( I ), Y( 1, I ), 1 )
274: *
275: * Compute T(1:i,i)
276: *
277: CALL ZSCAL( I-1, -TAU( I ), T( 1, I ), 1 )
278: CALL ZTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T, LDT,
279: $ T( 1, I ), 1 )
280: T( I, I ) = TAU( I )
281: *
282: 10 CONTINUE
283: A( K+NB, NB ) = EI
284: *
285: RETURN
286: *
287: * End of ZLAHRD
288: *
289: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlahrd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack