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