Annotation of rpl/lapack/lapack/zlahrd.f, revision 1.15
1.12 bertrand 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.
1.9 bertrand 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: *>
1.15 ! bertrand 37: *> This routine is deprecated and has been replaced by routine ZLAHR2.
! 38: *>
1.9 bertrand 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: *
1.15 ! bertrand 125: *> \date November 2015
1.9 bertrand 126: *
127: *> \ingroup complex16OTHERauxiliary
128: *
129: *> \par Further Details:
130: * =====================
131: *>
132: *> \verbatim
133: *>
134: *> The matrix Q is represented as a product of nb elementary reflectors
135: *>
136: *> Q = H(1) H(2) . . . H(nb).
137: *>
138: *> Each H(i) has the form
139: *>
140: *> H(i) = I - tau * v * v**H
141: *>
142: *> where tau is a complex scalar, and v is a complex vector with
143: *> v(1:i+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in
144: *> A(i+k+1:n,i), and tau in TAU(i).
145: *>
146: *> The elements of the vectors v together form the (n-k+1)-by-nb matrix
147: *> V which is needed, with T and Y, to apply the transformation to the
148: *> unreduced part of the matrix, using an update of the form:
149: *> A := (I - V*T*V**H) * (A - Y*V**H).
150: *>
151: *> The contents of A on exit are illustrated by the following example
152: *> with n = 7, k = 3 and nb = 2:
153: *>
154: *> ( a h a a a )
155: *> ( a h a a a )
156: *> ( a h a a a )
157: *> ( h h a a a )
158: *> ( v1 h a a a )
159: *> ( v1 v2 a a a )
160: *> ( v1 v2 a a a )
161: *>
162: *> where a denotes an element of the original matrix A, h denotes a
163: *> modified element of the upper Hessenberg matrix H, and vi denotes an
164: *> element of the vector defining H(i).
165: *> \endverbatim
166: *>
167: * =====================================================================
1.1 bertrand 168: SUBROUTINE ZLAHRD( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY )
169: *
1.15 ! bertrand 170: * -- LAPACK auxiliary routine (version 3.6.0) --
1.1 bertrand 171: * -- LAPACK is a software package provided by Univ. of Tennessee, --
172: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.15 ! bertrand 173: * November 2015
1.1 bertrand 174: *
175: * .. Scalar Arguments ..
176: INTEGER K, LDA, LDT, LDY, N, NB
177: * ..
178: * .. Array Arguments ..
179: COMPLEX*16 A( LDA, * ), T( LDT, NB ), TAU( NB ),
180: $ Y( LDY, NB )
181: * ..
182: *
183: * =====================================================================
184: *
185: * .. Parameters ..
186: COMPLEX*16 ZERO, ONE
187: PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ),
188: $ ONE = ( 1.0D+0, 0.0D+0 ) )
189: * ..
190: * .. Local Scalars ..
191: INTEGER I
192: COMPLEX*16 EI
193: * ..
194: * .. External Subroutines ..
195: EXTERNAL ZAXPY, ZCOPY, ZGEMV, ZLACGV, ZLARFG, ZSCAL,
196: $ ZTRMV
197: * ..
198: * .. Intrinsic Functions ..
199: INTRINSIC MIN
200: * ..
201: * .. Executable Statements ..
202: *
203: * Quick return if possible
204: *
205: IF( N.LE.1 )
206: $ RETURN
207: *
208: DO 10 I = 1, NB
209: IF( I.GT.1 ) THEN
210: *
211: * Update A(1:n,i)
212: *
1.8 bertrand 213: * Compute i-th column of A - Y * V**H
1.1 bertrand 214: *
215: CALL ZLACGV( I-1, A( K+I-1, 1 ), LDA )
216: CALL ZGEMV( 'No transpose', N, I-1, -ONE, Y, LDY,
217: $ A( K+I-1, 1 ), LDA, ONE, A( 1, I ), 1 )
218: CALL ZLACGV( I-1, A( K+I-1, 1 ), LDA )
219: *
1.8 bertrand 220: * Apply I - V * T**H * V**H to this column (call it b) from the
1.1 bertrand 221: * left, using the last column of T as workspace
222: *
223: * Let V = ( V1 ) and b = ( b1 ) (first I-1 rows)
224: * ( V2 ) ( b2 )
225: *
226: * where V1 is unit lower triangular
227: *
1.8 bertrand 228: * w := V1**H * b1
1.1 bertrand 229: *
230: CALL ZCOPY( I-1, A( K+1, I ), 1, T( 1, NB ), 1 )
231: CALL ZTRMV( 'Lower', 'Conjugate transpose', 'Unit', I-1,
232: $ A( K+1, 1 ), LDA, T( 1, NB ), 1 )
233: *
1.8 bertrand 234: * w := w + V2**H *b2
1.1 bertrand 235: *
236: CALL ZGEMV( 'Conjugate transpose', N-K-I+1, I-1, ONE,
237: $ A( K+I, 1 ), LDA, A( K+I, I ), 1, ONE,
238: $ T( 1, NB ), 1 )
239: *
1.8 bertrand 240: * w := T**H *w
1.1 bertrand 241: *
242: CALL ZTRMV( 'Upper', 'Conjugate transpose', 'Non-unit', I-1,
243: $ T, LDT, T( 1, NB ), 1 )
244: *
245: * b2 := b2 - V2*w
246: *
247: CALL ZGEMV( 'No transpose', N-K-I+1, I-1, -ONE, A( K+I, 1 ),
248: $ LDA, T( 1, NB ), 1, ONE, A( K+I, I ), 1 )
249: *
250: * b1 := b1 - V1*w
251: *
252: CALL ZTRMV( 'Lower', 'No transpose', 'Unit', I-1,
253: $ A( K+1, 1 ), LDA, T( 1, NB ), 1 )
254: CALL ZAXPY( I-1, -ONE, T( 1, NB ), 1, A( K+1, I ), 1 )
255: *
256: A( K+I-1, I-1 ) = EI
257: END IF
258: *
259: * Generate the elementary reflector H(i) to annihilate
260: * A(k+i+1:n,i)
261: *
262: EI = A( K+I, I )
263: CALL ZLARFG( N-K-I+1, EI, A( MIN( K+I+1, N ), I ), 1,
264: $ TAU( I ) )
265: A( K+I, I ) = ONE
266: *
267: * Compute Y(1:n,i)
268: *
269: CALL ZGEMV( 'No transpose', N, N-K-I+1, ONE, A( 1, I+1 ), LDA,
270: $ A( K+I, I ), 1, ZERO, Y( 1, I ), 1 )
271: CALL ZGEMV( 'Conjugate transpose', N-K-I+1, I-1, ONE,
272: $ A( K+I, 1 ), LDA, A( K+I, I ), 1, ZERO, T( 1, I ),
273: $ 1 )
274: CALL ZGEMV( 'No transpose', N, I-1, -ONE, Y, LDY, T( 1, I ), 1,
275: $ ONE, Y( 1, I ), 1 )
276: CALL ZSCAL( N, TAU( I ), Y( 1, I ), 1 )
277: *
278: * Compute T(1:i,i)
279: *
280: CALL ZSCAL( I-1, -TAU( I ), T( 1, I ), 1 )
281: CALL ZTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T, LDT,
282: $ T( 1, I ), 1 )
283: T( I, I ) = TAU( I )
284: *
285: 10 CONTINUE
286: A( K+NB, NB ) = EI
287: *
288: RETURN
289: *
290: * End of ZLAHRD
291: *
292: 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