LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
clatrd.f
Go to the documentation of this file.
1*> \brief \b CLATRD reduces the first nb rows and columns of a symmetric/Hermitian matrix A to real tridiagonal form by an unitary similarity transformation.
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download CLATRD + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clatrd.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clatrd.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clatrd.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE CLATRD( UPLO, N, NB, A, LDA, E, TAU, W, LDW )
20*
21* .. Scalar Arguments ..
22* CHARACTER UPLO
23* INTEGER LDA, LDW, N, NB
24* ..
25* .. Array Arguments ..
26* REAL E( * )
27* COMPLEX A( LDA, * ), TAU( * ), W( LDW, * )
28* ..
29*
30*
31*> \par Purpose:
32* =============
33*>
34*> \verbatim
35*>
36*> CLATRD reduces NB rows and columns of a complex Hermitian matrix A to
37*> Hermitian tridiagonal form by a unitary similarity
38*> transformation Q**H * A * Q, and returns the matrices V and W which are
39*> needed to apply the transformation to the unreduced part of A.
40*>
41*> If UPLO = 'U', CLATRD reduces the last NB rows and columns of a
42*> matrix, of which the upper triangle is supplied;
43*> if UPLO = 'L', CLATRD reduces the first NB rows and columns of a
44*> matrix, of which the lower triangle is supplied.
45*>
46*> This is an auxiliary routine called by CHETRD.
47*> \endverbatim
48*
49* Arguments:
50* ==========
51*
52*> \param[in] UPLO
53*> \verbatim
54*> UPLO is CHARACTER*1
55*> Specifies whether the upper or lower triangular part of the
56*> Hermitian matrix A is stored:
57*> = 'U': Upper triangular
58*> = 'L': Lower triangular
59*> \endverbatim
60*>
61*> \param[in] N
62*> \verbatim
63*> N is INTEGER
64*> The order of the matrix A.
65*> \endverbatim
66*>
67*> \param[in] NB
68*> \verbatim
69*> NB is INTEGER
70*> The number of rows and columns to be reduced.
71*> \endverbatim
72*>
73*> \param[in,out] A
74*> \verbatim
75*> A is COMPLEX array, dimension (LDA,N)
76*> On entry, the Hermitian matrix A. If UPLO = 'U', the leading
77*> n-by-n upper triangular part of A contains the upper
78*> triangular part of the matrix A, and the strictly lower
79*> triangular part of A is not referenced. If UPLO = 'L', the
80*> leading n-by-n lower triangular part of A contains the lower
81*> triangular part of the matrix A, and the strictly upper
82*> triangular part of A is not referenced.
83*> On exit:
84*> if UPLO = 'U', the last NB columns have been reduced to
85*> tridiagonal form, with the diagonal elements overwriting
86*> the diagonal elements of A; the elements above the diagonal
87*> with the array TAU, represent the unitary matrix Q as a
88*> product of elementary reflectors;
89*> if UPLO = 'L', the first NB columns have been reduced to
90*> tridiagonal form, with the diagonal elements overwriting
91*> the diagonal elements of A; the elements below the diagonal
92*> with the array TAU, represent the unitary matrix Q as a
93*> product of elementary reflectors.
94*> See Further Details.
95*> \endverbatim
96*>
97*> \param[in] LDA
98*> \verbatim
99*> LDA is INTEGER
100*> The leading dimension of the array A. LDA >= max(1,N).
101*> \endverbatim
102*>
103*> \param[out] E
104*> \verbatim
105*> E is REAL array, dimension (N-1)
106*> If UPLO = 'U', E(n-nb:n-1) contains the superdiagonal
107*> elements of the last NB columns of the reduced matrix;
108*> if UPLO = 'L', E(1:nb) contains the subdiagonal elements of
109*> the first NB columns of the reduced matrix.
110*> \endverbatim
111*>
112*> \param[out] TAU
113*> \verbatim
114*> TAU is COMPLEX array, dimension (N-1)
115*> The scalar factors of the elementary reflectors, stored in
116*> TAU(n-nb:n-1) if UPLO = 'U', and in TAU(1:nb) if UPLO = 'L'.
117*> See Further Details.
118*> \endverbatim
119*>
120*> \param[out] W
121*> \verbatim
122*> W is COMPLEX array, dimension (LDW,NB)
123*> The n-by-nb matrix W required to update the unreduced part
124*> of A.
125*> \endverbatim
126*>
127*> \param[in] LDW
128*> \verbatim
129*> LDW is INTEGER
130*> The leading dimension of the array W. LDW >= max(1,N).
131*> \endverbatim
132*
133* Authors:
134* ========
135*
136*> \author Univ. of Tennessee
137*> \author Univ. of California Berkeley
138*> \author Univ. of Colorado Denver
139*> \author NAG Ltd.
140*
141*> \ingroup latrd
142*
143*> \par Further Details:
144* =====================
145*>
146*> \verbatim
147*>
148*> If UPLO = 'U', the matrix Q is represented as a product of elementary
149*> reflectors
150*>
151*> Q = H(n) H(n-1) . . . H(n-nb+1).
152*>
153*> Each H(i) has the form
154*>
155*> H(i) = I - tau * v * v**H
156*>
157*> where tau is a complex scalar, and v is a complex vector with
158*> v(i:n) = 0 and v(i-1) = 1; v(1:i-1) is stored on exit in A(1:i-1,i),
159*> and tau in TAU(i-1).
160*>
161*> If UPLO = 'L', the matrix Q is represented as a product of elementary
162*> reflectors
163*>
164*> Q = H(1) H(2) . . . H(nb).
165*>
166*> Each H(i) has the form
167*>
168*> H(i) = I - tau * v * v**H
169*>
170*> where tau is a complex scalar, and v is a complex vector with
171*> v(1:i) = 0 and v(i+1) = 1; v(i+1:n) is stored on exit in A(i+1:n,i),
172*> and tau in TAU(i).
173*>
174*> The elements of the vectors v together form the n-by-nb matrix V
175*> which is needed, with W, to apply the transformation to the unreduced
176*> part of the matrix, using a Hermitian rank-2k update of the form:
177*> A := A - V*W**H - W*V**H.
178*>
179*> The contents of A on exit are illustrated by the following examples
180*> with n = 5 and nb = 2:
181*>
182*> if UPLO = 'U': if UPLO = 'L':
183*>
184*> ( a a a v4 v5 ) ( d )
185*> ( a a v4 v5 ) ( 1 d )
186*> ( a 1 v5 ) ( v1 1 a )
187*> ( d 1 ) ( v1 v2 a a )
188*> ( d ) ( v1 v2 a a a )
189*>
190*> where d denotes a diagonal element of the reduced matrix, a denotes
191*> an element of the original matrix that is unchanged, and vi denotes
192*> an element of the vector defining H(i).
193*> \endverbatim
194*>
195* =====================================================================
196 SUBROUTINE clatrd( UPLO, N, NB, A, LDA, E, TAU, W, LDW )
197*
198* -- LAPACK auxiliary routine --
199* -- LAPACK is a software package provided by Univ. of Tennessee, --
200* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
201*
202* .. Scalar Arguments ..
203 CHARACTER UPLO
204 INTEGER LDA, LDW, N, NB
205* ..
206* .. Array Arguments ..
207 REAL E( * )
208 COMPLEX A( LDA, * ), TAU( * ), W( LDW, * )
209* ..
210*
211* =====================================================================
212*
213* .. Parameters ..
214 COMPLEX ZERO, ONE, HALF
215 parameter( zero = ( 0.0e+0, 0.0e+0 ),
216 $ one = ( 1.0e+0, 0.0e+0 ),
217 $ half = ( 0.5e+0, 0.0e+0 ) )
218* ..
219* .. Local Scalars ..
220 INTEGER I, IW
221 COMPLEX ALPHA
222* ..
223* .. External Subroutines ..
224 EXTERNAL caxpy, cgemv, chemv, clacgv, clarfg,
225 $ cscal
226* ..
227* .. External Functions ..
228 LOGICAL LSAME
229 COMPLEX CDOTC
230 EXTERNAL lsame, cdotc
231* ..
232* .. Intrinsic Functions ..
233 INTRINSIC min, real
234* ..
235* .. Executable Statements ..
236*
237* Quick return if possible
238*
239 IF( n.LE.0 )
240 $ RETURN
241*
242 IF( lsame( uplo, 'U' ) ) THEN
243*
244* Reduce last NB columns of upper triangle
245*
246 DO 10 i = n, n - nb + 1, -1
247 iw = i - n + nb
248 IF( i.LT.n ) THEN
249*
250* Update A(1:i,i)
251*
252 a( i, i ) = real( a( i, i ) )
253 CALL clacgv( n-i, w( i, iw+1 ), ldw )
254 CALL cgemv( 'No transpose', i, n-i, -one, a( 1, i+1 ),
255 $ lda, w( i, iw+1 ), ldw, one, a( 1, i ), 1 )
256 CALL clacgv( n-i, w( i, iw+1 ), ldw )
257 CALL clacgv( n-i, a( i, i+1 ), lda )
258 CALL cgemv( 'No transpose', i, n-i, -one, w( 1,
259 $ iw+1 ),
260 $ ldw, a( i, i+1 ), lda, one, a( 1, i ), 1 )
261 CALL clacgv( n-i, a( i, i+1 ), lda )
262 a( i, i ) = real( a( i, i ) )
263 END IF
264 IF( i.GT.1 ) THEN
265*
266* Generate elementary reflector H(i) to annihilate
267* A(1:i-2,i)
268*
269 alpha = a( i-1, i )
270 CALL clarfg( i-1, alpha, a( 1, i ), 1, tau( i-1 ) )
271 e( i-1 ) = real( alpha )
272 a( i-1, i ) = one
273*
274* Compute W(1:i-1,i)
275*
276 CALL chemv( 'Upper', i-1, one, a, lda, a( 1, i ), 1,
277 $ zero, w( 1, iw ), 1 )
278 IF( i.LT.n ) THEN
279 CALL cgemv( 'Conjugate transpose', i-1, n-i, one,
280 $ w( 1, iw+1 ), ldw, a( 1, i ), 1, zero,
281 $ w( i+1, iw ), 1 )
282 CALL cgemv( 'No transpose', i-1, n-i, -one,
283 $ a( 1, i+1 ), lda, w( i+1, iw ), 1, one,
284 $ w( 1, iw ), 1 )
285 CALL cgemv( 'Conjugate transpose', i-1, n-i, one,
286 $ a( 1, i+1 ), lda, a( 1, i ), 1, zero,
287 $ w( i+1, iw ), 1 )
288 CALL cgemv( 'No transpose', i-1, n-i, -one,
289 $ w( 1, iw+1 ), ldw, w( i+1, iw ), 1, one,
290 $ w( 1, iw ), 1 )
291 END IF
292 CALL cscal( i-1, tau( i-1 ), w( 1, iw ), 1 )
293 alpha = -half*tau( i-1 )*cdotc( i-1, w( 1, iw ), 1,
294 $ a( 1, i ), 1 )
295 CALL caxpy( i-1, alpha, a( 1, i ), 1, w( 1, iw ), 1 )
296 END IF
297*
298 10 CONTINUE
299 ELSE
300*
301* Reduce first NB columns of lower triangle
302*
303 DO 20 i = 1, nb
304*
305* Update A(i:n,i)
306*
307 a( i, i ) = real( a( i, i ) )
308 CALL clacgv( i-1, w( i, 1 ), ldw )
309 CALL cgemv( 'No transpose', n-i+1, i-1, -one, a( i, 1 ),
310 $ lda, w( i, 1 ), ldw, one, a( i, i ), 1 )
311 CALL clacgv( i-1, w( i, 1 ), ldw )
312 CALL clacgv( i-1, a( i, 1 ), lda )
313 CALL cgemv( 'No transpose', n-i+1, i-1, -one, w( i, 1 ),
314 $ ldw, a( i, 1 ), lda, one, a( i, i ), 1 )
315 CALL clacgv( i-1, a( i, 1 ), lda )
316 a( i, i ) = real( a( i, i ) )
317 IF( i.LT.n ) THEN
318*
319* Generate elementary reflector H(i) to annihilate
320* A(i+2:n,i)
321*
322 alpha = a( i+1, i )
323 CALL clarfg( n-i, alpha, a( min( i+2, n ), i ), 1,
324 $ tau( i ) )
325 e( i ) = real( alpha )
326 a( i+1, i ) = one
327*
328* Compute W(i+1:n,i)
329*
330 CALL chemv( 'Lower', n-i, one, a( i+1, i+1 ), lda,
331 $ a( i+1, i ), 1, zero, w( i+1, i ), 1 )
332 CALL cgemv( 'Conjugate transpose', n-i, i-1, one,
333 $ w( i+1, 1 ), ldw, a( i+1, i ), 1, zero,
334 $ w( 1, i ), 1 )
335 CALL cgemv( 'No transpose', n-i, i-1, -one, a( i+1,
336 $ 1 ),
337 $ lda, w( 1, i ), 1, one, w( i+1, i ), 1 )
338 CALL cgemv( 'Conjugate transpose', n-i, i-1, one,
339 $ a( i+1, 1 ), lda, a( i+1, i ), 1, zero,
340 $ w( 1, i ), 1 )
341 CALL cgemv( 'No transpose', n-i, i-1, -one, w( i+1,
342 $ 1 ),
343 $ ldw, w( 1, i ), 1, one, w( i+1, i ), 1 )
344 CALL cscal( n-i, tau( i ), w( i+1, i ), 1 )
345 alpha = -half*tau( i )*cdotc( n-i, w( i+1, i ), 1,
346 $ a( i+1, i ), 1 )
347 CALL caxpy( n-i, alpha, a( i+1, i ), 1, w( i+1, i ),
348 $ 1 )
349 END IF
350*
351 20 CONTINUE
352 END IF
353*
354 RETURN
355*
356* End of CLATRD
357*
358 END
caxpy
subroutine caxpy(n, ca, cx, incx, cy, incy)
CAXPY
Definition caxpy.f:88
cgemv
subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
CGEMV
Definition cgemv.f:160
chemv
subroutine chemv(uplo, n, alpha, a, lda, x, incx, beta, y, incy)
CHEMV
Definition chemv.f:154
clacgv
subroutine clacgv(n, x, incx)
CLACGV conjugates a complex vector.
Definition clacgv.f:72
clarfg
subroutine clarfg(n, alpha, x, incx, tau)
CLARFG generates an elementary reflector (Householder matrix).
Definition clarfg.f:104
clatrd
subroutine clatrd(uplo, n, nb, a, lda, e, tau, w, ldw)
CLATRD reduces the first nb rows and columns of a symmetric/Hermitian matrix A to real tridiagonal fo...
Definition clatrd.f:197
cscal
subroutine cscal(n, ca, cx, incx)
CSCAL
Definition cscal.f:78